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VARIATIONAL  ANALYSIS  AND  APPROXIMATE 
SOLUTIONS  FOR  TRANSPORT  PHENOMENA 

INTRODUCTION 

Approximate  methods  have  been  applied  to  the  solution  of  partial  differential  equations  when 
analytical  forms  of  solution  are  not  possible  due  to  the  nature  of  the  particular  problems  under  study. 
Some  types  of  approximate  solutions  can  be  directly  applied  to  the  given  equation  and  some  require 
first  the  derivation  of  a  variational  form  of  this  equation.  The  latter  imposes  certain  restrictions  on 
applying  approximate  methods  since  many  of  the  variational  derivations  are  problem  dependent. 

The  intent  of  this  study  is  to  provide  the  formulation  of  a  unified  analysis  for  partial  differential 
equations,  such  as  the  general  transport  equation,  by  methods  analogous  to  those  of  classical  mechan¬ 
ics.  The  implementation  of  certain  concepts  of  classical  mechanics  opens  the  way  to  a  formulation  of 
the  transport  equation  by  means  of  generalized  coordinates  and  leads  to  Lagrangian  type  of  equations. 
This  approach  simplifies  the  application  of  approximate  methods  for  obtaining  solutions  to  many  practi¬ 
cal  problems  of  interest. 

The  variational  analysis  presented  here  involves  a  generalization  of  the  concept  of  virtual  work 
from  classical  mechanics.  The  approach  was  first  adapted  by  M.A.  Biot  and  was  applied  to  heat  transfer 
phenomena  and  irreversible  processes  [1,2].  By  extending  the  concept  of  virtual  work  to  the  general 
transport  equation  a  variational  form  of  this  equation  can  be  obtained,  which  is  problem  independent, 
and  the  resulting  variational  equation  contains  terms  similar  to  those  found  in  Lagrangian  mechanics. 
This  approach  should  not  be  regarded  as  a  "new  method"  for  deriving  variational  formulations.  It  is 
only  an  application  of  classical  concepts  and  an  extension  of  Lagrangian  mechanics  to  derive  a  unified 
formulation  for  the  transport  equation.  One  of  the  advantages  of  such  approach  is  that  the  physical 
significance  of  each  term  of  the  differential  equation  can  be  easily  identified. 
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The  governing  equation  for  transport  phenomena  and  the  introduction  of  a  vector  defined  as  tran¬ 
sport  displacement  are  given  in  the  first  section  of  this  study.  Based  on  the  definition  of  the  displace¬ 
ment  vector  the  transport  equation  is  written  in  a  form  similar  to  the  equation  of  motion  in  mechanics. 
In  this  form  the  governing  equation  will  be  used  in  the  next  section  to  derive  the  variational  formula¬ 
tion  for  the  transport  equation. 

According  to  concepts  of  classical  mechanics  two  types  of  variational  formulations  can  be  derived, 
the  fundamental  form  based  on  variations  of  the  displacement  field  and  the  complementary  form  based 
on  variations  of  the  stress  field.  The  basic  definitions  and  derivations  for  the  variational  analysis  are 
given  in  the  second  section  and  the  two  types  of  variational  formulations  are  presented  for  the  transport 
equation.  The  analogy  between  the  derived  equations  and  those  of  mechanics  can  be  justified  by  identi¬ 
fying  each  of  the  terms  of  the  variational  equations  to  be  similar  to  such  quantities  as  potential  func¬ 
tion,  dissipation  function  and  generalized  body  and  boundary  forces.  The  dissipation  function  is  a  gen¬ 
eralization  of  the  concept  introduced  by  Rayleigh  in  mechanics  for  systems  with  viscous  dissipation. 
The  resulting  Lagrangian  equations,  in  terms  of  generalized  coordinates,  represented  either  type  of  vari¬ 
ational  formulation  and  they  are  not.  problem  dependent.  As  a  consequence  the  generalized  equations 
can  be  used  for  obtaining  solutions  to  physical  problems  governed  by  such  equations  as  the  transport 
equation.  Furthermore,  many  types  of  approximate  solutions  can  be  derived  by  expressing  the  field 
parameters  in  terms  of  the  generalized  coordinates. 

A  particular  example  is  given  in  the  third  section  where  a  series  expansion  is  assumed  for  the  field 
parameters.  This  expansion  expresses  the  field  parameters  as  a  linear  combination  of  the  generalized 
coordinates.  Although  such  an  approximation  is  restricted  by  the  linear  dependence  on  the  generalized 
coordinates,  it  may  represent  such  approximations  as  a  Fourier  series  expansion  or  an  expansion  in 
terms  of  orthogonal  functions.  The  derivation  of  the  equations  in  terms  of  the  above  approximation  is 
given  and  it  is  shown  that  the  finite  element  methods  can  be  derived  as  a  special  case  of  this  approxi¬ 
mation  [31.  The  discretized  system  of  equations  for  a  linear  finite  element  approximation  is  given  as  an 
example  and  four  types  of  finite  element  models  are  derived.  Two  of  the  models  correspond  to  the 
fundamental  variational  formulation  and  the  other  two  correspond  to  the  complementary  one. 
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Numerical  experiments  and  solutions  to  particular  problems  are  given  in  the  fourth  section  of  this 
study.  In  the  first  part  of  the  numerical  applications  the  error  behavior  of  the  finite  element  models  is 
investigated  for  a  given  problem  with  known  analytical  solutions.  Furthermore,  the  convergence  of  the 
models  is  investigated  and  the  criteria  for  uniform  convergence  are  discussed.  In  the  second  part  of  the 
numerical  applications,  examples  of  numerical  solutions  are  presented  for  a  number  of  physical  prob¬ 
lems.  From  the  example  presented  the  efficiency  of  the  computational  models  is  discussed  and  one  can 
observe  that  the  unified  formulation  presented  here  can  be  extended  to  a  variety  of  physical  problems. 


1.  BASIC  EQUATIONS  AND  DEFINITIONS 


The  general  form  of  a  partial  differential  equation  describing  transport  phenomena,  such  as  mass 
transfer,  heat  transfer  or  species  concentration  is  given  by 


dC(xk,t) 


dx, 


+  S  -  K  (C  -  C0) 


(1) 


where  C  is  the  transport  variable,  x,  is  the  coordinate  system,  t  is  the  time,  V,  is  the  transport  velocity 
in  x,  direction,  ktJ  is  the  tensor  dispersion  coefficient,  S  is  the  generalized  source  term  and  K  is  the  gen¬ 
eralized  transfer  coefficient. 


In  order  to  describe  the  kinematics  of  a  medium,  reference  is  made  to  two  configurations.  The 
reference  or  initial  configuration  at  time  t  —  t„  and  the  present  configuration  at  time  t  with  C„  and  C 
the  values  of  the  transport  variable  respectively.  A  dimensionless  transport  variable  can  be  defined  in 
terms  of  C  and  C0  if  one  considers  the  ratio  of  the  change  AC  of  the  absolute  value  of  the  transport 
variable  from  the  reference  to  the  present  state  over  the  absolute  value  of  that  variable  at  the  reference 
state 


@(xk.t) 


C-C„ 


AC 

c„  • 


In  terms  of  the  dimensionless  variable  0,  Eq.  (1)  is  expressed  as 


(2) 


(3) 
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From  the  physical  point  of  view  the  dimensionless  variable  0  can  be  interpreted  as  a  generalized 
deformation  similar  to  mechanical  strain.  Define  now  a  vector  field  H,  (xk,t)  as  the  transport  displace¬ 
ment  vector ■,  which  is  a  function  of  the  coordinates  xk  and  time  r,  of  the  form 


and 


H,(xk,t) 


dff,(xk,t) 

Ft 


(4) 


&(xk,t) 


dHj(xk,t) 

a.v, 


(5) 


The  interpretation  of  0  as  a  deformation  or  dilatation  (strain)  is  justified  by  Eq.  (5).  Further¬ 
more,  Eq.  (5)  may  be  considered  as  a  constraint  in  the  sense  of  classical  mechanics  and  it  must  be 
verified  by  the  physical  solution  and  by  the  variations  of  8 H,  and  80  as 

8©  - (3//,).  (6) 

OXj 

The  variations  SH,  correspond  to  virtual  displacements  with  0  and  //,  being  analogous  to  strain  and  dis¬ 
placement  in  mechanics.  Another  quantity,  which  is  defined  in  mechanics  and  it  is  related  to  deforma¬ 
tion,  is  the  stress  tensor.  Here  a  similar  stress  can  be  defined  as 


cr  -  £0  (7) 

where  £  is  a  modulus  similar  to  the  buik  modulus  in  mechanics.  Eq.  (7)  may  represent  a  constitutive 

relation  for  a  particular  medium  and  together  with  Eq.  (3)  describes  the  behavior  of  that  medium  under 

certain  loading  conditions. 


The  governing  equation  can  be  expressed  in  terms  of  the  displacement  field  as 


dt 


where  satisfied  the  equation 


d  H,  30 

+  KQ-  ku  ff-  -  KH,  +  h, 


(8) 


6h,  1  ft 

-T--  TrJn  Sdt.  (9) 

dx,  C0  J  o 

This  equation  does  not  uniquely  determine  h,  but  any  particular  field  may  be  chosen  to  satisfy  Eq.  (9). 
Such  a  field  may  be  considered  to  be  a  given  function  of  the  coordinates  and  time.  The  dispersion 
coefficient  is  a  tensor  with  six  components  and  with  the  property 
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k,j  “  kj! , 

hence,  ku  is  a  symmetric  tensor.  If  Ay  is  the  inverse  of  ku  then  A,,  is  also  symmetric  and  Eq.  f 8)  can 
be  written  as 


BH, 


Bt 


+  K© 


A© 

+  k,J(KH,-hl)  -  f~ 


(10) 


Equations  (5)  and  (8)  represent  a  formulation  which  is  defined  as  the  displacement  formulation  for  the 
transport  equation  and  this  formulation  reduces  to  the  conventional  one,  Eq.  (3),  by  combining  Eqs. 
(5)  and  (8).  Furthermore,  the  three  equations,  Eqs.  (5),  (7)  and  (8),  are  analogous  to  the  kinematic 
relations,  stress-strain  relations  and  the  momentum  equation  in  mechanics. 


The  concept  of  transport  displacement  has  been  introduced  previously  by  the  author,  under  the 
name  of  heat  displacement,  for  describing  heat  transfer  processes  [41.  In  general,  transport  processes 
can  be  described  either  by  the  conventional  transport  equation,  Eq.  (3),  or  by  the  two  equivalent  equa¬ 
tions.  Eqs.  (5)  and  (8).  The  introduced  transport  displacement  is  regarded  as  a  generalized  quantity 
conjugate  to  the  generalized  deformation  ©.  This  generalized  presentation  of  the  governing  equation 
for  transport  phenomena  has  certain  features  which  will  be  discussed  in  the  next  sections. 

2.  VARIATIONAL  ANALYSIS* 

A  large  number  of  variational  formulations  have  been  introduced  in  the  past  for  deriving  approxi¬ 
mate  solutions.  The  majority  of  them  are  based  on  the  minimization  of  a  functional  which  describes  a 
particular  physical  process.  Such  formulations  are  usually  restricted  to  the  conditions  of  the  particular 
problem  and  their  applicability  is  limited.  One  wishes  to  have  a  unified  approach  for  deriving  varia¬ 
tional  formulations.  Such  formulations  should  be  problem  independent  and  they  should  be  derived 
from  physical  considerations. 

Classical  mechanics  has  provided  researchers  with  some  powerful  tools  for  solving  problems  either 
analytically  or  by  approximate  methods.  One  such  a  tool  is  the  principle  of  virtual  work  which  has  been 
successful  in  solving  a  variety  of  problems.  The  principle  of  virtual  work  in  mechanics  has  been 
extended  to  such  areas  as  thermodynamics  [5]  and  variational  equations  can  be  derived  which  are  not 
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problem  dependent  for  a  variety  of  physical  processes.  Another  concept  largely  used  in  mechanics  is 
the  concept  of  generalized  coordinates  which  can  be  used  to  describe  a  physical  system.  Generalized 
coordinates  can  be  interpreted  as  generalized  displacements  or  velocities  which  specify  the 
configurations  of  a  physical  system  [6], 

The  application  of  these  concepts  are  not  restricted  to  mechanics  or  only  to  certain  media.  In  this 
section  the  concept  of  virtual  work  is  used  to  derive  two  variational  formulations  for  the  transport  equa¬ 
tion.  The  first  is  referred  to  as  the  fundamental  form  and  is  based  on  the  displacement  formulation  of 
the  transport  equation.  The  second  formulation  is  based  on  the  conventional  form  of  the  transport 
equation  and  is  referred  to  as  the  complementary  form.  These  two  forms  are  directly  analogous  to  the 
ones  in  classical  mechanics  where  the  fundamental  form  of  a  variational  principle  is  given  in  terms  of 
variations  of  the  displacement  field  and  the  complementary  form  in  terms  of  forces.  Both  of  these 
forms  are  based  on  equilibrium  or  conservation  laws  governing  the  physical  system. 

2.1  Fundamental  Variational  Principle 

The  fundamental  form  is  derived  from  the  displacement  formulation  of  the  transport  equation  by 
implementing  the  principle  of  virtual  work.  Consider  a  variation  8 Ht  of  the  displacement  field  H,  and 
the  corresponding  variation  80  given  by  the  constraint  Eq.  (6).  If  the  medium  is  subjected  to  a  virtual 
displacement  8//,,  then  the  principle  of  virtual  work  requires  that 

-t^+V,©  +\JKH,  -  h,)  -  6H,ctv-Q.  (11) 

v  or  oXj 

where  the  integration  is  extended  over  a  volume  v  of  the  medium.  Integration  of  Eq.  (11)  by  parts 
yields 

j>„  K&  8//,</v+j' ’©  -^~^H,)dv+J K.jiKH,- hjfSH.dv  =  f  Qn:bH,ds.  (12) 

The  surface  integral  is  extended  over  the  boundary  surface  s  of  the  control  volume  v  and  n,  is  the  unit 
vector  normal  to  the  boundary  pointing  outward.  The  second  term  of  the  left  hand  side  is  written  as 
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(13) 


r  oh.  r 

I  08  — —  i/v  =  I  080c/v  =  8P 

•'v  O.V,  J  v 

and  the  scalar  function  P  may  be  interpreted  as  a  potential  function,  similar  to  the  potential  function  of 
mechanics,  given  by 

P  =  1/2  f  92dv.  (14) 

In  terms  of  the  variations  8 P,  Eq.  (12)  is  written  as 

/<  A  //  i*  /• 

Jv  ~  +  Vj9  8Hidv  +  Jv  xw(^  -  V  8 H t dv  +  8P  -  J  9n,8H,  ds.  (15) 

The  variational  relation  given  by  Eq.  (15)  can  be  considered  as  the  variational  principle  for  the  tran¬ 
sport  equation.  Equation  (15)  must  be  verified  for  arbitrary  variations  of  the  displacement  field  H,. 
with  0  defined  as  a  function  of  H )  through  the  constraint,  Eq.  (5).  Hence  the  variational  equation 
given  by  Eq.  (15)  is  nothing  but  a  variational  statement  for  the  transport  equation  with  energy  conser¬ 
vation  automatically  satisfied.  This  variational  principle  must  be  understood  in  a  broad  sense  since  it 
corresponds  to  the  principle  of  virtual  work  in  mechanics.  It  represents  the  fundamental  form  since  it 
was  derived  in  terms  of  the  displacement  field. 


2.2  Complementary  Variational  Principle 

In  the  preceding  section  the  fundamental  form  of  the  variational  principle  was  derived  for  the 
transport  equation.  The  analogy  with  mechanics  suggests  that  it  is  possible  to  derive  a  complementary 
form  of  the  variational  principle  for  the  same  equation.  Since  the  generalized  deformation  0,  Eq.  (5). 
leads  to  the  definition  of  the  generalized  stress  cr,  which  is  similar  to  the  mechanical  stress,  then  the 
complementary  form  is  derived  by  varying  the  stress  field  and  expressing  the  conservation  equation, 
Eq.  (3),  in  variational  form. 


In  terms  of  the  generalized  stress  o-,  Eq.  (3)  is  written  as 


IT  +  T~  ( v‘a)  ~  T~  k'i  T~  +  Ka  ~  s  =  0 

dt  dx,  ox,  dXj  C0E 

and  for  arbitrary  variation  8<r  an  equivalent  to  Eq.  (11)  is  obtained  as 


+  Ka-  - 


8c rc/v  =*  0. 


(16) 


(17) 


Integration  by  parts  yields 


f 

8cr</v+f  k a  (ba)dv  -t-  f 

A  aba  S5<j 

Jv 

dt  dx,  ‘ 

Jw  1  dx,  dx, 

CnE 

d\ 


C  /  9<t  j 

k,,n,  — —  as. 
1  1  dx. 


(18) 


The  surface  integral  is  extended  to  the  boundary  surface  s  of  the  volume  v  and  n,  denotes  the  unit  vec¬ 
tor  normal  to  the  boundary  pointing  outwards.  Proceeding  as  in  the  previous  section  Eq.  (18)  is  written 


as 


f 

iz.  a.  \ 

badv-t  f 

Ka  -  —5—  5 

badv  +  bP  =  f 

J  V 

dr 

'  dx, 

J  V 

C0E  j 

J  $ 

3tr 


(19) 


where  the  variation  bP  is  given  by 


f 

ktJ¥~b 

da 

dv, 

-) — - —  aba 

J  V 

"  dXj 

dx, 

dx, 

d\ 


(20) 


and  the  scalar  P  bv 


>«/. 


SK 


,  da  da _ . 

,J  dx ,  dx,  dx,  <J 


d\. 


(21) 


Equation  (19)  represents  the  complementary  form  of  the  variational  principle  for  the  transport  equa¬ 
tion.  The  different  terms  of  the  variational  equation,  Eq.  (19)  may  be  evaluated  either  in  terms  of  the 
generalized  stress  a  or  in  terms  of  the  deformations  ®.  The  latter  case  corresponds  to  formulating  the 
variational  principle  for  the  transport  equation  in  terms  of  a  primitive  dimensionless  variable  0. 


The  two  forms  of  variational  principles  derived  previously  are  generalized  statements  of  conserva¬ 
tion  laws  for  a  physical  system  and  their  existence  is  verified  by  concepts  from  classical  mechanics. 
Either  form  can  be  applied  to  the  governing  equation  of  many  physical  processes,  but  it  first  is  neces¬ 
sary  to  derive  a  more  compact  form  of  these  variational  equations  which  will  be  better  suited  to  practi¬ 
cal  applications.  Some  remarks  can  be  made  regarding  the  form  of  the  variational  equations.  A 
significant  difference  between  the  two  forms  is  the  fact  that  the  complementary  form  involves  space 
derivatives  of  the  transport  variable  0.  This  does  not  present  a  problem  as  long  as  the  space  distribu¬ 
tion  for  0  is  continuous.  If  discontinuities  exist,  as  is  the  case  in  many  practical  problems,  then  the 
complementary  form  is  not  the  appropriate  one  to  be  used  for  the  solution  of  such  a  problem.  In  gen- 
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eral.  when  using  the  same  approximate  representation  for  the  transport  variable  0,  the  complementary 
principle  will  yield  results  which  are  less  accurate  than  those  obtained  by  the  fundamental  form  which  is 
based  on  the  displacement  field. 

2.3  Generalized  Coordinates 

A  physical  system  is  usually  described  by  a  set  of  coordinates  and  the  governing  equations  are 
expressed  in  terms  of  these  coordinates.  In  many  cases,  a  given  set  of  coordinates  is  not  the  most  con¬ 
venient  one  to  describe  a  physical  problem.  As  an  example,  for  certain  problems  involving  the  motion 
of  a  particle  it  is  more  convenient  to  use  a  moving  coordinate  system  than  a  fixed  one.  It  would  be 
very  desirable,  and  practical,  to  have  a  general  method  for  describing  a  physical  system  by  a  set  of  coor¬ 
dinates  which  will  provide  a  uniform  method  to  represent  and  solve  the  governing  equations  for  physi¬ 
cal  processes.  Such  a  set  is  given  by  the  generalized  coordinates  of  classical  mechanics.  A  set  of  gen¬ 
eralized  coordinates  is  any  set  of  coordinates  by  means  of  which  the  position  of  a  particle  in  a  system 
may  be  specified.  The  concept  of  the  generalized  coordinates  and  the  method  for  describing  a  physical 
system  in  terms  of  these  coordinates  was  first  introduced  by  Lagrange. 

When  a  physical  system  is  described  by  a  set  of  generalized  coordinates,  it  is  the  usual  practice  to 

designate  each  coordinate  by  q  and  the  set  of  n  generalized  coordinates  bv  qk,  k  =  1 . n.  For  a 

physical  system  of  .V  particles  a  set  of  3.V  generalized  coordinates  may  be  specified,  and  since  these 
coordinates  must  have  some  definite  set  of  values  they  will  be  functions  of  the  time  and  also  ol  the 
cartesian  system  of  coordinates.  For  a  physical  system  which  is  described  by  a  set  of  generalized  coordi¬ 
nates  qk,  the  time  derivative  qk  is  defined  as  a  generalized  velocity  associated  with  this  coordinate.  This 
generalized  velocity,  for  example,  can  be  computed  in  terms  of  cartesian  coordinates  and  velocities.  If 
one  assumes  that  qk  is  a  function  of  the  cartesian  coordinates  x,  and  time  r.  then 

Qk  =  Qk f  -v, .  Z  > 

It  is  also  possible  to  express  the  cartesian  coordinates  in  terms  of  generalized  coordinates  as 

-v,  -  x,{qk.t) 

From  the  definition  of  the  generalized  velocity  one  obtains 
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w 


X, 


5-v,  dx,  . 

"s7  +  Qk 


Other  physical  quantities,  such  as  energy,  momentum,  forces  etc,  can  be  expressed  in  terms  of 
the  generalized  coordinates.  For  example  the  generalized  momentum  pk  associated  with  the  coordinate 
qk  is  given  by 

d  T 

Pk  =  — 

OQk 

when  Tis  the  kinetic  energy  in  terms  of  qk. 


The  subject  of  this  section  is  to  implement  the  concept  of  generalized  coordinates  to  the  two 
forms  of  variational  equations  previously  derived.  By  introducing  generalized  coordinates  into  the  vari¬ 
ational  equations,  they  can  be  translated  into  a  Lagrangian  type  of  equations.  This  provides  a  unified 
approach  and  the  transport  equation  is  represented  in  a  generalized  form  independent  of  coordinate  sys¬ 
tems  or  physical  conditions. 


Consider  the  displacement  field  H,  to  be  a  function  of  a  set  of  n  generalized  coordinates  qk< 

H,  -  H,  (qk,  x„  ,).  (22) 

The  variations  of  Hn  or  small  displacements  8 H,,  are  due  entirely  to  the  variations  of  the  generalized 

coordinates  dqk  and  they  are  expressed  as 


dH, 

oH,  =>  — —  8qk. 
oqk 


(23) 


The  variations  5qk  can  be  identified  as  virtual  displacements  since  they  do  not  necessarily  represent  any 
actual  motion  but  only  any  possible  motion  of  a  system.  Under  such  a  motion  there  is  an  amount  of 
work  done  which  is  defined  as  the  virtual  work  8  IT  and  it  is  expressed  as 


««'-0*8fc  (24) 

where  Qk  is  defined  as  a  generalized  force  associated  with  the  coordinate  qK. 


The  variations  of  the  potential  function  P,  given  by  Eq.  (14),  are  expressed  as 


5P  “  8^ 

oqk 

and  the  time  derivative  of  the  displacement  vector  is  taken  as 


(25) 
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126) 


H  - 


dH, 

dqk 


<h  + 


dH, 

at 


Equation  (26)  expresses  the  total  time  derivative  of  the  vector  //,  and  by  differentiating  with 
respect  to  the  generalized  velocity  qk  one  obtains 


an,  sh_ 
aqk  3^ 


(27) 


Taking  into  account  Eqs.  (23),  (25),  (26)  and  (27),  the  variational  principle,  Eq.  (15).  can  be 
transformed  into  a  more  compact  form.  The  first  term  of  Eq.  (15)  is  written  as 


f 

dH, 

- -  +  f,0 

dH.dv  =»  f  A . . 

IJ 

dt 

/  1  1/ 

dt 

3W,  ,  .  dD 

— —  a qkdv  -  — —  6(7* 
dqk  dqk 


where  the  function  D  is  given  by 


D  -  1/2  /  \„(H,  +  (,<-))  (/V,  +  r,(-))c/v  -  1/2  f  K„H’H’dv 
and  the  vector  H*  is  given  by 


(28) 


(29) 


//,*=//,+  F,0.  (30) 

The  vector  H'  can  be  regarded  as  the  total  rate  of  the  transport  displacement  and  the  function  D 

defines  a  dissipation  function.  A  physical  interpretation  of  both  of  these  qualities  can  be  given  if  one 
considers  some  particular  case  of  a  physical  process.  For  example,  if  the  transport  equation  describes 
heat  transfer,  then  the  vector  H‘  represents  the  total  convective  and  diffusive  local  rate  of  heat  flow 
and  the  function  D  the  total  dissipation  of  the  physical  process.  This  dissipation  has  been  previously 
introduced  for  systems  with  friction  or  viscous  effects  and  it  is  known  as  Rayleigh's  dissipation  func¬ 
tion. 


In  terms  of  the  variations  8 P  and  8 D.  Eq.  (15)  yields 


8  D 


C  3//, 

+  dP  =*  J  Wtj,  — —  8 qk  ds  + 
Js  oqk 


8F. 


(31) 


The  term  8F of  the  right  hand  side  of  Eq.  (31)  combines  the  source  and  transfer  terms  and  is  given  by 

d  H 

8F-  f  \J-KH,  +  ih)  8</*,/v.  (32) 

•'v  "  rtj. 
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and 


H,H,  +  hjH, 


dv 


(33) 


The  surface  integral  of  the  right  hand  side  of  Eq.  (31)  can  be  expressed  in  terms  of  a  generalized  boun¬ 
dary  force  Qk  as 


r  on, 

Js  -7^-  &Qkds  “  Q^dk 

where 


(34) 


r  on, 

Qk  “  J  &V,  t — -  rfs.  (35) 

•/s 

From  the  physical  point  of  view,  Eq.  (34),  can  be  interpreted  as  the  virtual  work,  SIT,  done  by  the 
force  Q,  on  the  inflow  of  a  volume  of  fluid  per  unit  area  associated  with  the  virtual  displacement  6 qk 
Comparing  Eqs.  (24),  (31)  and  (34)  the  relation 


8D  4- SP  -  Sf  =  SIT  (36) 

is  obtained.  This  can  be  regarded  as  expressing  conservation  since  the  virtual  work  done  by  the  gen¬ 
eralized  boundary  forces  is  equal  to  energy  gain  or  loss  for  a  physical  system. 


For  arbitrary  variations  of  the  generalized  coordinates  Eq.  (36)  is  written  as 

w+lL-Qk  +  JLL. 

dqk  37  k  37* 


(37) 


For  the  set  of  n  generalized  coordinates  Eqs.  (37)  constitute  a  system  of  n  differential  equation  for  the 
unknowns  qk.  These  equations  are  of  the  same  form  as  the  equations  of  Lagrangian  mechanics  for  the 
slow  motion  of  a  dissipative  system  with  negligible  inertia  forces.  For  example,  for  the  particular  case 
of  heat  transfer  the  potential  function  is  regarded  as  a  thermal  potential  and  the  function  D  as  the  total 
heat  dissipation  function.  The  function  Pmay  be  due  to  internal  heat  generation  and/or  heat  transfer. 
Furthermore,  the  generalized  forces  Qk  represent  thermal  forces  due  to  the  temperature  distribution  at 
the  boundary  and  they  are  defined  by  Eq.  (34)  as  the  work  done  by  Q,  under  the  displacement  8qk. 


This  generalized  derivation  in  terms  of  the  generalized  coordinates,  can  be  extended  to  the  com¬ 
plementary  form  of  the  variational  principle.  The  derivation  is  identical  and  it  yields  a  set  of  equations 
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similar  to  Eqs.  (37),  the  difference  being  that  they  are  derived  in  terms  of  the  transport  variable  0.  If 
one  assumes  that  0  is  a  function  of  a  set  of  n  generalized  coordinates  qk,  then 


0  -  0(v*.  -V„  /).  (38) 

Following  the  same  steps  as  for  the  previous  derivations  the  Lagrangian  equations  are 


where 


M  + J'.fc  +  iZ. 

Bq  Bqk  Bqk 


(39) 


D 

P 


112 1, 
1/2 1, 


0  +  V 


,,  90_ 
'  Bx, 


0  +  V, 


50. 

Bx, 


dv 


,  50  50  d  K  , 


dv 


-1/2  K 02  +  ~~  SQ 


d\ 


(40) 

Ul) 

(42) 


and 


<?*  “  X  k»  9 


50  50 
5.v,  dqk 


ds. 


(43) 


As  one  may  observe  from  Eqs.  (40)  to  (43)  the  complementary  variational  principle  translates  to  simi¬ 
lar  types  of  equations  as  the  fundamental  form,  but  there  is  an  important  difference  between  the  tv  o 
formulations.  The  latter  form  involves  spatial  derivatives  of  the  transport  variable  0  in  the  expressions 
for  the  functions  D ,  P  and  the  generalized  forces  Q,.  This  is  a  disadvantage  of  the  complementary 
form,  as  it  was  mentioned  earlier,  especially  for  applications  to  problems  with  discontinuities  on  the  dis¬ 
tribution  of  the  primitive  variables.  Although  each  of  the  formulations  was  derived  in  terms  of  a 
different  set  of  coordinates  both  can  be  represented  by  the  same  type  of  generalized  equations.  This  is 
verified  by  the  unified  approach  used  to  derive  the  Lagrangian  type  of  equations  and  also  by  the  fact 
that  the  transport  variable  0  and  the  displacement  field  are  conjugate  variables  related  by  the  constraint 
Eq.  (5). 


The  concept  of  generalized  coordinates  provides  the  basis  for  deriving  the  Lagrangian  formulation 
of  the  transport  equation.  From  the  physical  point  of  view  the  generalized  coordinates  can  describe  a 
system  completely  and  the  Lagrangian  equations  can  provide  an  accurate  formulation  of  the  physical 


13 


problem.  A  particular  advantage  of  the  formulation  is  that  it  was  derived  independently  of  a  particular 
representation  of  the  field  variables.  Furthermore,  the  above  analysis  can  be  extended  to  many  physical 
systems  governed  by  equations  other  than  the  transport  equation. 

The  Lagrangian  types  of  equations  derived  in  this  section  do  not  represent  a  new  mathematical 
theory  but  a  different  way  of  expressing  a  conservation  law.  The  main  feature  of  these  equations  is  the 
manner  in  which  they  were  derived.  It  is  evident  that  these  equations  have  the  same  form  in  any  sys¬ 
tem  of  coordinates  and  from  classical  mechanics  one  can  verify  that  the  functions  involved  in 
Lagrange’s  equations  have  the  same  value  for  any  coordinate  system.  The  unified  form  of  these  equa¬ 
tions  is  of  great  value  not  only  for  advanced  formulations  in  mechanics  but  they  also  provide  a  starting 
point  for  approximate  solutions  to  physical  problems. 

3.  APPROXIMATE  METHODS 

The  partial  differential  equation  describing  transport  phenomena  was  transformed  into  a  set  of 
ordinary  differential  equations  expressed  in  terms  of  generalized  coordinates.  The  equations  obtained 
can  be  solved  analytically  for  a  number  of  physical  problems.  For  many  of  the  problems  encountered 
in  practical  applications,  analytical  solutions  are  not  easy  to  obtain  and  approximate  or  numerical  solu¬ 
tions  are  sought. 

The  generalized  nature  of  the  previously  obtained  equations  is  most  appropriate  for  applying 
approximate  solutions.  Since  these  equations  are  given  in  terms  of  generalized  coordinates,  the  field 
variables  describing  the  physical  process  are  assumed  to  be  functions  of  these  generalized  coordinates. 
Many  different  expressions  or  combinations  of  the  generalized  coordinates  can  be  used  to  form  the 
functional  representing  the  field  variable.  In  order  to  demonstrate  the  procedure  one  may  assume  a 
linear  dependence  of  the  field  variables  on  the  generalized  coordinates.  This  specific  representation  by 
no  means  restricts  the  previously  derived  formulations  to  only  linear  functional  approximation.  The 
main  reason  for  such  a  choice  is  the  fact  that  a  variety  of  numerical  solutions  are  based  in  such  approxi¬ 
mations  and  as  will  be  shown  later  the  finite  element  method  can  be  derived  as  a  special  case  of  the 
generalized  formulations. 
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3.1  Linear  Functionals 


The  linear  dependence  of  the  displacement  field  H,  on  the  generalized  coordinates  is  expressed  as 

H,(qk,  Xj.  t)  =  fik(Xj)qk(t)  (44) 

The  generalized  coordinates  qk  may  be  functions  of  time  and  they  are  treated  as  the  coordinates 

describing  the  physical  system.  From  the  physical  point  of  view  they  represent  degrees  of  freedom  and 
the  functions  fk  (.v,)  specify  the  extent  to  which  qk  participate  in  the  functions  H,Lxj.t).  The  approxi¬ 
mation  of  the  displacement  field  given  by  Eq.  (44)  is  quite  general  since  no  specific  expressions  are 
given  for  the  functions  flk.  Furthermore,  Eq.  (44)  may  represent  approximations  such  as  Fourier  series 
or  series  expansions  in  terms  of  orthogonal  functions.  Similar  approximations  can  be  used  for  the  tran¬ 
sport  variable  0. 


In  the  following  the  fundamental  form  of  the  variational  principle  will  be  used  to  derive  the 
approximate  system  of  equations  for  obtaining  numerical  solutions.  For  the  complementary  form  the 
derivation  will  not  be  given  since  it  is  identical  to  that  of  the  fundamental  form. 


The  Lagrangian  type  of  equation  given  by  Eq.  (37)  can  be  translated  to  a  system  of  ordinary 
differential  equations  in  terms  of  the  generalized  coordinates  qk  by  using  the  approximation,  Eq.  (44), 
for  the  displacement  field  and  by  first  defining  the  time  and  spatial  derivatives  as 

ffi  ~  Qkfki  f45) 

and 


on, 

®  -  it  • 

In  terms  of  these  definitions  the  dissipation  function  D  is  evaluated  as 

D  “  1/2  l  ^.Mkfki  +  K %/*«.*) <<?////  +  VjqiJ)nn)dv 
or 

D  -  1/2  d,Jq,q)  +  \,Jq,qJ 

where 


(46) 


d'J  “  X  *  klJkifl/dV 


(47) 
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and 


v-/“J/w^Wv'  <48> 

Similarly,  the  potential  function  Z3  is  evaluated  as 

P  ~  \/2  p,Jq,q)  (49) 

where 

P„  =*  fv  fk.kfju^-  (50) 

The  function  F  transforms  to  the  form 


F  =  S,q,  -  \/2wlt  q,  q; 

where 

S,-f  XWA,. fads  (51) 

and 

wq  “  /v  “  S^ufikfjidv.  (52) 

Finally  the  boundary  forces  Q,  are  expressed  as 

Qi  “  fs&^j/ijds  (53) 

Substituting  the  expressions  for  D,  Py  F,  and  Q,  from  the  above  equations  into  Eq.  (44)  yields 

d,,Qj  +  (v,,  +  p,t  +  w^qj  =  Q,.  +  S,.  (54) 

This  matrix  equation  represents  a  system  or  ordinary  differential  equations  for  the  unknown  field 
parameters  qh  which  may  represent  displacements.  The  generalized  coordinates  describe  the  physical 
system  by  means  of  a  large  but  finite  number  of  variables  and  their  number  will  determine  the  size  of 
the  above  system  of  equations.  The  solution  of  Eq.  (54)  yields  the  approximate  solution  for  the  tran¬ 
sport  equation.  The  order  of  approximation  of  the  numerical  solution  is  determined  by  the  functions 
f,i.\k)  and  depends  on  the  accuracy  of  the  solution  required  and  usually  on  the  complexity  of  the  phy¬ 
sical  problem. 

For  the  evaluation  of  the  matrix  coefficients  of  Eq.  (54)  the  integration  domain  is  considered  to 
be  taken  over  the  volume  v  in  the  rt-space  of  the  coordinates  .v,,  where 
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d\  =  dx  i  dx 2  ...  dx„. 

In  a  more  general  sense  the  integration  domain  can  be  extended  over  all  n  +  1  coordinates  .v,  and  t 
instead  of  only  the  n  coordinates  .v,,  where 


dr  =  dx\,  dxi  . . .  dxn  dt. 

This  extension  of  the  integration  domain  will  yield  the  same  Lagrangian-type  equations  since  the  unk¬ 
nown  physical  quantities  can  be  expressed  in  any  functional  form  in  terms  of  generalized  coordinates. 
For  this  particular  case  the  approximation  given  by  Eq.  (44)  will  be  replaced  by 


H,(xk.t)  -  4  <**> &(/)«;.  (55) 

It  should  be  noted  that  the  approximation  given  by  Eq.  (55)  is  more  restrictive  in  nature  than  the  one 

given  by  Eq.  (44)  since  one  assumes  a  known  dependence  of  the  variables  on  time  through  the  func¬ 
tion  g„(t).  The  generalized  coordinates,  in  this  case  q".  represent  values  at  the  i  th  point  in  space  and 
at  the  n  th  point  in  time.  The  respective  time  and  spatial  derivatives  of  Eq.  (55)  yield 


H,  -  fij(xk)  gn(t)  q; 
hu  -  f,jj  (**)  sn  (')?;. 


(56) 


A  system  of  equations  in  terms  of  the  generalized  coordinates  q"  can  be  obtained  by  substituting  the 
approximation  given  by  Eq.  (55)  into  the  Lagrangian  equations.  The  expressions  for  D,  P ,  Fand  Q, 
are  evaluated  as  for  the  previous  approximations  and  they  are 


D  =  1/2  ldijk,  +  v„J  qfqj 

P  -  \llplMqW)  (57) 

F  =  S,kq^  -  1/2 iv,. ^,4' 

where 

dtjkl  ~  ^  mnJ tmf  jn&k 

„  V ijkl  *  § r  ^  mn  ^ mJ tnfjt.sKk  Rl 

Pijkl  -f,  f.m.mfjn.ngkSldT  (58) 

^:k  “  k mn^nfrmgk^7 

wijkl  *  ^  mnftmf jngkgld7  ■ 

Substituting  the  above  expressions  into  the  Lagrangian  Eq.  (44)  yields 


K /*/  +  V./U  +  Piikl  +  ^jkiUt1 


Q!  +  S„ 


i.  j  “  1  ,n 
k.l  —  1  ,m 


(59) 
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where  the  boundary  forces  Qj  are  given  by 

Qj  ~  fs  Vi/km.m/,jgl8k  (60) 

The  matrix  Eq.  (59)  represents  a  system  of  n  +  m  algebraic  equations  for  the  unknown  field  parame¬ 
ters,  the  generalized  coordinates  qj,  and  the  domain  of  integration  r  of  the  matrix  coefficients  is 
extended  over  the  n  +  1  coordinates  x,  and  t. 

The  Lagrangian  type  of  equations  and  the  two  systems  of  differential  and  algebraic  equations 
derived  in  this  section  represent  two  examples  of  approximate  methods  which  can  be  applied  to  solve 
problems  on  transport  processes.  The  application  of  the  derived  formulation  to  a  physical  system  can 
be  carried  out  by  dividing  the  domain  of  the  system  into  a  finite  number  of  subdomains  connected  by 
common  boundaries.  Each  of  these  domains  constitutes  a  sub-system  that  can  be  analyzed  separately 
by  the  system  of  equations  derived  in  terms  of  generalized  coordinates.  For  each  sub-system  s  the  two 
systems  of  equations  may  be  written  in  a  general  form  as 

4s)q,  -  Q>s)>  (61) 

where  q,  are  the  generalized  coordinates  for  the  sub-system  s,  Q;(s>  are  the  associated  generalized  boun¬ 
dary  forces  and  Ajjs)  the  matrix  coefficients  of  the  system  of  equations  which  were  derived  from  the 
Lagrangian  equation  for  the  domain  s.  Taking  the  sum  of  Eq.  (61)  over  all  the  domains  of  the  sub¬ 
systems,  the  equation  for  the  total  system  is  written  as 

Avqj  -  Q„  (62) 

where  Ay  is  now  the  matrix  coefficient  for  the  system  of  equations  corresponding  to  the  total  domain 

and  Qj  are  the  forces  only  at  the  boundary  of  the  total  system.  This  result  is  justified  by  the  fact  that 
the  summation  extends  to  the  variables  not  only  at  the  interfaces  of  the  sub-systems  but  also  to  the 
variables  at  the  boundary  of  the  total  system.  However,  at  the  interfaces  the  forces  are  in  pairs  that 
cancel  out.  For  example,  for  two  sub-systems  s  and  s' with  a  common  boundary  the  generalized  coordi¬ 
nates  at  the  interface  are  q,  and  the  corresponding  forces  are  Q,  and  Q,.  From  the  definition  of  the 
boundary  force,  the  outward  normal  vectors  for  the  domains  with  common  boundaries  are  in  opposite 
directions,  hence 

q,  -  -  q: 
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or 


(63) 


Q,  +  Qi  -  o. 

This  leads  to  the  cancellation  of  all  the  forces  at  the  interfaces  and  the  derived  Eqs.  (62)  include  only 
the  forces  at  the  boundary  of  the  total  system. 

The  procedure  of  dividing  the  domain  of  a  system  to  sub-spaces  is  commonly  used  in  approximate 
methods  and  a  typical  example  is  the  identical  procedure  used  for  the  analysis  of  a  physical  system  by 
the  finite  element  method. 

3.2  Finite  Element  Method 

A  variety  of  approaches  can  be  found  under  the  name  of  the  finite  element  method.  Such 
approaches  include  the  well  known  Galerkin  methods,  the  Rayleigh-Ritz  methods  and  many  others 
which  are  usually  modifications  or  derivatives  of  the  ones  mentioned.  The  formulation  and  derivation 
of  the  Galerkin  method  are  well  known  and  the  method  has  been  successfully  applied  to  a  variety  of 
problems  for  obtaining  numerical  solutions.  The  Galerkin  approach  is  usually  referred  to  as  the  con¬ 
ventional  finite  element  method  in  contrast  to  those  which  are  based  on  modification  of  the  Galerkin  or 
non-conventional  approaches.  The  intention  here  is  to  show  that  the  conventional  finite  element 
method  is  a  special  application  of  the  present  generalized  formulation.  Consider  the  following  approxi¬ 
mation  for  the  displacement  field  H, 

H,(x,t)  —  a0  +  <7|.v  +  a2x 2  +  a}x}  +  . . .  .  (64) 

and  the  corresponding  deformation  field  © 

0  =  at  +  2a2x  +  3a3.v:  +  . . .  .  (65) 

where  the  a,  U  —  0,  1,  . . .  rt)  are  time  dependent  coefficients  to  be  determined.  Although  the  approxi¬ 
mation  given  by  Eq.  (64)  is  restricted  to  one  dimensional  space,  it  can  be  easily  extended  to  two  or 
three  dimensions  since  the  derived  formulations  for  approximate  methods  are  of  general  form.  In 
order  to  demonstrate  the  derivation  of  the  matrix  equation  for  the  finite  element  method  only  a  first 
order  approximation  will  be  considered  here.  Higher  order  approximations  have  been  considered  in 
previous  studies  [3,4]  and  their  derivation  is  carried  out  in  a  similar  way  as  the  first  order  one.  For  a 


I  9 


first  order  approximation  the  derived  finite  element  models  correspond  to  linear  element  models  with 
two  degrees  of  freedom.  Four  such  linear  elements  will  be  considered  here.  The  first  two  are  linear 
displacement  models  denoted  by  DFEM  and  DFET  and  correspond  for  the  approximations  given  by 
Eqs.  (44)  and  (55).  The  last  two  are  conventional  linear  models  in  terms  of  the  primitive  variable  0 
and  are  denoted  by  CFEM  and  CFET.  The  displacement  models  are  derived  from  the  fundamental 
form  of  the  variational  formulation  and  the  conventional  models  from  the  complementary  form. 

3.2.1  Displacement  Models 

(a)  DFEM 


From  Eq.  (64)  the  approximation  of  the  displacement  field  Hl.x.t)  for  the  linear  element  of 
length  A.v  is  given  by 


HLx.t) 


1  - 


A.v 


Hx(t)  +  H2U) 


(66) 


where  H\  and  H2  are  the  nodal  values  of  the  displacement.  They  can  be  identified  as  the  generalized 
coordinates  qx  and  q2  with  the  corresponding  basis  functions  fxl  and  fx2  given  by 


f,  i 


1  - 


A.v 


and  J'x2 


A.v 


Then,  according  to  Eq.  (44),  the  approximation  Eq.  (66)  is 


HLx.t)  =  fx,(x) 

For  the  deformation  0  the  approximation  is  evaluated  from  Eq.  (67)  as 

0U,r)  -  Wt)  -  Hx(t)). 

A.v 


(67) 


(68) 


(69) 


Note  that  within  each  element  0  varies  only  with  time  for  the  linear  approximation  of  the  displacement 
field.  This  is  analogous  to  a  constant  strain  element  in  mechanics. 


The  matrix  coefficients  of  Eq.  (54)  are  evaluated  in  terms  of  Eqs.  (66)  and  (69)  as 


[v„ 


.4  A.v 
6  k 


2  1 

,  ,  4 

1  -l 

1  2 

■ 

-1  1 

-1  1 
■1  1 


(70) 


K,] 


2  1 
I  2 


A  V 
2  k 
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AK\.x 
6  k 


and 


Q,  I 


q i  =  a  -  &  -  .40 

t-a.v 


where  A  is  the  cross  sectional  area  of  the  element,  k  is  the  diffusion  coefficient  and  V  the  transport 
velocity.  If  one  also  assumes  a  linear  approximation  for  the  source  term  then  S,  from  Eq.  (51)  is 
evaluated  as 


A  A.v 
6  k 


2  1 

ki 

A  A.v 

f-s. 

1  2 

6  k 

l* 

(71) 


where  5)  =  2h\  +  fu  and  Si  =  h\  +  2h->.  The  matrix  equation  for  the  DFEM  model  is  then  assembled 
as 


2  1 

Hi 

3FU 

-1  l 

6 

1  -1 

2  1 

N 

k  ! 

[-©] 

[5. 

1  2 

Hi 

+ 

K 

-I  1 

+  K 

-1  1 

+  K0 

1  2 

N 

fl 

1  ©1 

4- 

^2, 

(72) 

The  matrix  equation,  Eq.  (72),  has  been  derived  in  dimensionless  form  according  to  the  transforma¬ 
tions 

C„  i 

-,-lh 

(73) 


x 

X~T 


L- 


T'.  © 


C  -  C  - 
— — -f,  //  = 

Ci  —  Ci 


W  -  ^  V  K 

o  L  ■  0  k  ■  *o 


KL 2  =  L 

,  t  ~  t . 


Cl  -  c0  L 

c„ 


k  Cx-  C 


Here  L  is  a  characteristic  length  and  C,  a  constant  reference  value  for  C.  The  bar  has  been  eliminated 
from  Eq.  (72)  for  simplicity. 

(b)  DFET  model 


If  one  now  considers  the  approximation  given  by  Eq.  (55),  then  the  linear  approximation  for 
Hix.t)  and  0(.v,r)  will  be  given  as 


H(x.t) 


1  -  T“ 

1  - 

Hi  + 

1  -  — 

—  H{  +  -f- 

Ax 

Ar 

A.v 

Ar  A.v 

1  Ar 

H>  +  -7—  -7-  Hi  (74) 
A.v  Ar 
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The  two  sets  of  equations  given  by  Eqs.  (72)  and  (78)  represent  the  displacement  finite  element 
models  for  obtaining  numerical  solutions  of  the  transport  equation.  Both  models  were  derived  from  the 
fundamental  form  of  the  variational  formulation  as  examples  of  approximate  methods. 

3. 2. 2  Con  ventional  Models 

(a)  CFE.VI  Model 

The  approximation  for  the  transport  variable  0  for  the  linear  element  of  length  A.v  is  given  by 

0(.v.r)  =  1  -  ~r~  0,  +  —0,  (79) 

A  a  A.v  ‘ 

where  0]  and  02  are  the  nodal  values  of  the  transport  variable.  They  can  be  identified  as  the  general¬ 
ized  coordinates  q\  and  <71  with  the  corresponding  functions  /vt  and  fsl  given  by  Eq.  (67). 

The  matrix  coefficients  and  the  expressions  for  D,  P,  and  F,  and  Q \  are  evaluated  from  Eqs.  (90) 
through  (43)  and  the  resulting  matrix  equation  is  given  by 
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2  1 

3i;, 

-1  1 

b 

i-i 

2  1 

f«i 

A 

~(0i 

-  0,) 

fs, 

1  2 

+ 

Wa 

-1  1 

-i  i 

+  Kn 

1  2 

U 

(0; 

-0.) 

The  above  equation  has  resulted  from  the  complementary  form  of  the  variational  formulation  and  is 
identical  to  the  finite  element  equations  obtained  by  the  conventional  Galerkin  method.  Furthermore 
the  matrix  coefficients  are  identical  to  the  ones  given  by  Eqs.  (70). 

(b)  CFET  Model 


The  linear  approximation  for  the  transport  variable  0  is  given  by 


0(.v,r)  = 


i 

1-  -L 

0/  + 

1  -  JL 

—  0,:  +  — 

\--L 

A.v 

lt 

A.v 

Si  1  A  a- 

lt 

0-!  +  —  —  ©r 
Ax  Ar  * 


(81) 


As  for  the  previous  approximations  the  generalized  coordinates  can  be  identified  as 


<7/  =  0/.  i  -  1.  2.  7=1.2  (82) 

and  the  corresponding  functions  /v,  and  g„  are  given  by  Eqs.  (76)  and  (77).  The  matrix  coefficients  are 
evaluated  in  terms  of  Eq.  (81)  and  the  matrix  equation  for  the  CFET  model  can  be  derived  from  Eq. 
(78)  by  replacing  the  nodal  unknowns  H>  by  0/. 


The  four  finite  element  models  presented  in  this  section  can  be  grouped  into  two  types  of  models 
which  correspond  to  the  general  approximation  given  by  Eqs.  (44)  and  (55).  The  first  type  can  be 
represented  by  the  following  generalized  equation 


d„q,  +  <v„  +  pu  +  wtJ )  q j  =  Q,  +  S,  (83) 

where  the  matrix  coefficients  are  given  by  Eqs.  (70)  and  the  generalized  coordinates  are  identified  as 


DFEM  model  :  ?,  =  //,  /  —  1,  2 
CFEM  model  :  <7,-0,  (=1,2. 

For  the  DFEM  model  Eq.  (83)  should  be  used  together  with  the  constraint 


(84) 


The  second  type  of  models  can  be  represented  by 

U„ki  +  v„ki  +  A,*/  +  W/m/1  q‘,  =  Q +  S* 


(85) 


(86) 
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where  the  matrices  are  given  by  Eqs.  (78)  and  the  generalized  coordinates  are  identified  as 


DFET  model;  q{  -  H-  i  -  I.  2.  /  -  !.  2 


(S'7) 


CFET  model:  q{  =  H;  <  =  1.  2,  j  =  1,2 

and  tor  the  DFET  model  Eq.  (83)  should  be  used  together  with  the  constraint  given  by  Eq.  (75).  In 
the  following  table  the  four  models  are  presented  in  summary. 


Table  1 


r - 

Approximation 

Model 

Generalized 

Coordinates 

Type  of  Matrix 
Equations 

i  Only  for 

spatial  coordinates 

DFEM 

q  =*  H 

Ordinary 

Differential 

S-type 

CFEM 

<7,  =  <->. 

Eq.  (83) 

I  For  both  spatial 

DFET 

<//  =  H; 

Algebraic 

i  and  time  coordinates 

Eq.  (86) 

•  S-T-type 

CFET 

_ 

<//  = 

Both  types  of  models  have  been  used  to  solve  problems  of  transport  phenomena  and  some  exam¬ 
ples  will  be  given  later  in  this  study.  From  these  examples  a  realistic  comparison  between  models  and 
also  an  evaluation  of  the  models  is  possible. 


The  generalized  variational  formulations  and  the  derived  approximate  methods  presented  here 
represent  a  unified  approach  for  describing  transport  phenomena.  Some  of  the  advantages  of  this 
unified  approach  are  as  follows: 

•  The  variational  formulation  is  problem  independent. 

•  The  presentation  of  the  transport  equation  by  generalized  coordinates  simplifies  the  complexity  of 
the  equation. 

•  The  Lagrangian  type  of  equations  are  most  suitable  for  obtaining  approximate  solutions. 

•  The  finite  element  method  can  be  derived  directly  as  a  special  application  for  approximate  solutions. 
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•  The  derived  finite  element  models  can  be  represented  by  one  generalized  model. 


•  The  displacement  models  are  most  suitable  for  solving  problems  with  discontinuities. 

In  the  following  a  number  of  numerical  solutions  will  be  presented  for  specific  boundary  value 
problems  and  the  error  behavior  of  these  numerical  solutions  will  be  investigated. 

4.  NUMERICAL  APPLICATIONS 


The  transport  equation  is  considered  in  this  section  as  describing  convective  heat  transfer  With 
this  choice  of  a  transport  process,  one  can  identify  the  transport  parameter  <-)  to  represent  temperature 


defined  by 


where  T  and  T„  are  absolute  temperature  at  the  present  and  reference  states  respectively  and  the  tran¬ 


sport  displacement  will  be  then  called  the  heat  displacement  [3], 


Convective  heat  transfer  provides  a  good  test  case  for  numerical  solutions  since  a  variety  of  boun¬ 
dary  conditions  can  be  considered.  Before  numerical  solutions  to  specific  problems  are  presented  it  is 
appropriate  to  investigate  the  stability  and  the  error  behavior  of  the  four  element  models.  Since  it  is 
necessary  to  use  a  numerical  scheme  to  approximate  the  time  derivatives  of  the  S-type  models,  a  finite 
difference  scheme  is  chosen  for  that  purpose.  The  scheme  is  known  as  the  backward  finite  difference 
approximation.  It  is  an  implicit,  unconditionally  stable  scheme  and  the  general  approximation  for  the 
first  derivative  is  given  by  the  following  series  expansion  [7] 


=  —  V  u,v<r,_,)  +  0( At"-1)  (88) 

Jt  fy  - 

where  a,  are  constant  coefficients  and  their  values  are  given  for  up  to  a  fourth  order  approximation  in 
Table  2. 


Table  2 


Order  of  Approximation 

«o 

«i 

a  2 

‘>4 

CI> 

First,  n  —  l 

1.0 

-1.0 

0.0 

0.0 

0.0 

0 

Second,  n  =  2 

1.5 

-2.0 

0.5 

0.0 

0.0 

0 

Third,  n  “  3 

11./6. 

-3.0 

1.5 

-1./3. 

0.0 

0 

Fourth,  n  =  4 

25/12. 

-4,0 

3.0 

-.4/3. 

0.25 

0 
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The  error  analysis  of  the  numerical  solution  of  a  particular  problem  will  be  presented  in  the  following 
section  and  some  criteria  for  convergence  and  error  estimates  will  be  discussed. 

4.1  Problem  Definition 


One  representative  example  of  transport  phenomena  is  the  convection  and  diffusion  of  heat  by  a 
moving  fluid.  The  governing  equation  for  the  one  dimensional  case  of  convection-diffusion  is  given  by 

t  *  r-  J7-d-0-°  18,1 

where  (■)  is  the  dimensionless  temperature.  I'„  is  the  constant  flow  velocity  and  D„  a  constant 
diffusivity.  It  has  been  also  assumed  that  there  is  no  heat  sources  or  sink  involved  in  the  fluid.  The 
boundary  and  initial  conditions  are  as  follows 


l-X.v.O)  =  0.0,  0.0  <  a 

0(0. r)  =  1.0.  0.0  iS  /  i,  (90) 

0<L.r)  =  0.0.  0.0  ^  t  ^  t, 

where  L  is  a  finite  characteristic  length  approximating  infinity.  The  analytical  solution  to  the  above 
problem  is  given  by 


<-)(.v j)  =*  1/2  exp 


y  +  V„r 

TJW 


v  -  K,t 

20D7' 


The  solution  of  the  convection-diffusion  equation  in  the  semi-infinite  space  is  a  typical  example  for 
testing  numerical  solutions  and  studying  error  behavior,  its  numerical  solution  is  challenging,  especially 
for  small  values  of  the  diffusion  coefficient  (D„  «  F(),  where  the  discontinuities  in  the  temperature 


field  are  difficult  to  simulate  numerically. 


4.2  Error  Analysis 


Numerical  solutions  of  the  above  problem  are  obtained  by  the  four  finite  element  models  previ¬ 
ously  derived  and  these  numerical  solutions  are  compared  to  the  analytical  one,  Eq.  (91),  in  order  to 
investigate  the  error  behavior.  For  the  S-type  models  (DFEM  and  CFEM)  a  first  order  backward  finite 
difference  scheme  is  chosen  to  approximate  the  time  derivative.  This  is  an  appropriate  choice  since  the 
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order  of  accuracy  of  the  first  order  scheme  is  the  same  as  the  accuracy  of  the  linear  in  time  approxima¬ 
tion  of  the  S-T-type  models  (DFET  and  CFET>. 


The  error  of  the  numerical  solution  is  the  average  error  oxer  the  entire  length  /.  at  a  particular 
time  level  and  i’  is  evaluated  by  the  following  relation 

E„  =  ~  T  !W  -  <-> .1  (92) 

where  .Vis  the  number  of  elements  used  for  the  space  discritization  and  (-)  -  ti,j  is  the  absolute  local 
error.  For  illustrating  numerical  results,  the  characteristic  length  L  =  1  is  divided  into  a  uniform  mesh 
with  .V£  elements  and  .VP  =  ,V£  +  1  nodal  points  (degrees  of  freedom),  and  the  dimensionless  ele¬ 
ment  length  then  equals  to  1/,V£.  Three  sets  of  values  for  the  constants  f',  and  D0  will  be  used  in 
investigating  the  convergence  and  error  of  the  numerical  solutions.  These  values  are  shown  in  Table  3 
and  each  set  is  used  for  all  four  models. 


Table  3 


Case 

D„ 

Name 

I 

■ 

0.001 

Convection 

II 

■SB 

1.0 

Convection-Diffusion 

III 

1.0 

Diffusion 

For  each  of  these  values  four  sets  of  results  will  be  presented  for  all  four  models.  The  first  set  of 
results  is  obtained  for  a  constant  value  of  A.v  and  variable  A/  (Table  4)  and  the  second  set  is  obtained 
for  constant  A /  and  variable  A.v  (Table  5).  The  third  for  constant  values  of  the  ratio  Az/A.v  (Table  6) 
and  the  fourth  set  for  constant  values  of  the  ratio  A//A.v;  (Table  7). 


Table  4 


NE 

A.v 

A  / 

V/,  -  A //A.v 

.V/:  -  A//A.v: 

0.002 

0.04 

0.00025 

0.005 

0.10 

0.010 

0.20 

0.025 

0.50 

0.100 

2.00 

0.200 

4.00 

0.250 

5.00 

1 

0.400 

8.00 

0.500 

10.00 

0.0500 

1.000 

20.00 

Table  5 


A  r 

NE 

Ax 

A/|  =»  Ar/Ax 

M1  =  Ar/Ax2 

5 

1/5.0 

mjssm 

10 

1/10.0 

mMmmzu, 

20 

1/20.0 

0.100 

30 

1/30.0 

0.150 

0.005 

40 

1/40.0 

50 

1/50,0 

— 

60 

1/60.0 

80 

1/80.0 

100 

200 

Table  6 


Ax 

Ar 

— 

A/,  —  Ar/Ax 

A/2  =  Ar/Ax2 

10 

0.100000 

0.02500 

2.50 

20 

0.05000 

0.01250 

5.00 

40 

0.02500 

0.00625 

10.00 

50 

0.02000 

0.00500 

0.25 

12.50 

80 

0.01250 

0.00312 

20.00 

100 

0.01000 

0.00250 

25.00 

Table  7 


NE 

Ax 

. 

Ar 

A/,  =  Ar/Ax 

M2  -  Ar/Ax2 

0.1 

0.250 

I 

0.00625 

0.125 

El 

1 

0.00156 

0.063 

2.5 

50 

1 

0.00100 

0.050 

80 

0.00040 

0.031 

100 

■IlfltSws! 

0.00025 

0.025 

In  all  of  the  figures  the  average  error  is  presented  at  three  time  levels.  This  is  necessary  in  order 
to  investigate  the  error  behavior  throughout  the  time  integration  domain.  The  first  set.  Figs.  1-3,  the 
value  of  is  constant  (NE  -  20)  and  the  time  step  Ar  was  given  values  0.0001  <  Ar  ^  0.05.  The 
corresponding  values  for  Ar/Ax  and  Ar/Ax2  are  given  in  Table  4.  The  second  set  of  results  are  given 
in  Figs.  4-6  where  A r  is  constant  (At  —  0.005)  and  NE  was  given  values  5  ^  NE  <  200.  Table  5  con¬ 
tains  the  values  for  Ax  and  the  corresponding  values  of  Ar/Ax  and  Ar/Ax2. 
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4.2.1  Displacement  Models  (DFET,  DFEM) 

In  Figs,  la  through  6a  results  for  the  average  errors  E„  are  given  for  the  two  displacement  models 
and  for  the  three  cases  of  the  coefficients  V0  and  D„.  In  the  first  three  figures  (la,  2a,  3a),  the  error  £„ 
is  given  as  a  function  of  It  and  the  next  three  figures,  (4a,  5a,  6a),  £„  is  given  as  a  function  of  ,V£. 
The  error  behavior  in  all  six  figures  follows  a  similar  pattern.  As  A /  or  ,V£  increases  there  is  initially  a 
reduction  in  the  error  of  the  numerical  solutions,  and  a  minimim  error  is  reached.  With  further 
increases  of  Ar  or  NE  the  error  starts  to  increase  but  the  rate  of  change  is  smaller  than  before  the 
minimum  was  reached.  The  difference  between  the  values  of  the  error  for  the  two  models  is  very  small 
for  small  values  of  A /  but  it  becomes  noticeable  as  A /  increases  with  the  DFET  model  being  the  more 
accurate  one.  Similar  observations  can  be  made  for  the  error  behavior  with  respect  to  V£.  The  error 
behavior  in  all  six  figures  shows  a  nonuniform  convergence  as  Af  increases  or  Ax  decreases.  This  pat¬ 
tern  is  not  unusual  since  only  A  t  or  Ax  is  changed  while  Ax  or  A /  remains  constant  respectively.  Cer¬ 
tain  observations  can  be  made  from  the  results  presented  in  the  six  figures. 

•  There  is  an  optimum  combination  of  Arand  Ax  for  minimum  error  for  both  models. 

•  To  the  left  of  the  optimum  point  the  error  of  the  DFEM  model  is  smaller  than  the  error  for  the 
DFET  model.  The  order  of  accuracy  is  reversed  just  before  the  optimum  point  and  the  difference 
in  error  between  the  two  models  becomes  larger. 

•  For  small  values  of  the  coefficient  D0 ,  ( D0  «  F0),  there  are  larger  differences  in  error  between 
the  two  models  with  the  DFET  being  more  accurate. 

•  As  time  progresses  the  value  of  the  minimum  error  decreases  and  the  optimum  points  are 
translated  to  the  right 

•  For  large  values  of  Aror  NE  the  error  increases  monotonically  with  time. 

•  There  seems  to  be  an  area  where  the  error  curves  for  different  time  levels  intersect,  indicating  that 
a  constant  error  value  can  be  achieved  for  certain  combinations  of  A  t  and  Ax  throughout  the  time 
domain.  This  is  more  pronounced  in  Figs.  4a-6a. 
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Fig.  1  —  Average  error  for  constant  Aar.  ( Y0  -  1.0.  D0  -  0.001) 
(a)  Displacement  models  and  (b)  Conventional  models 


30 


AVEK.  EkRQH 


*  CFET 


.  53L  -  CrtM 


Z.  00  .  3  '.  .32  .33  .  04 


T : vp  ~T£P  QX 

(hi 

Fig.  2  —  Average  error  for  constant  Ax  ( V0  -  1.0.  D0  —  1.0): 
(a)  Displacement  models  and  (b)  Conventional  models 
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Fig.  3  —  Average  error  for  constam  A*.  (  E0  -  0.0.  £?0  -  1.0): 
(a)  Displacement  models  and  (b)  Conventional  models 
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Fig.  4  —  Average  error  for  constant  Ar.  (  F0  -  1.0.  £>0  -  0.001): 

(a)  Displacement  models  and  (b)  Conventional  models 
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Fig.  5  —  Average  error  for  constant  Ar,  (  Va  -  1.0,  O0  -  1 
(a)  Displacement  models  and  (b)  Conventional  models 
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Fig.  6  —  Average  error  for  conslant  .ir.  ( F0  -  0.0,  D0  -  1 
(a)  Displacement  models  and  (b)  Conventional  models 
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The  similarities  observed  in  error  behavior  may  lead  to  some  optimum  combination  of  Ar  and  Ax  for 
error  control  of  the  numerical  solutions. 

4.2.2  Conventional  Models  (CFET.  CFEM) 

In  Figs.  lb-6b  results  are  given  for  the  two  conventional  finite  element  models  and  for  the  same 
set  of  values  of  the  coefficients  D0  and  V0  as  before.  The  error  behavior  with  respect  to  Ar  or  Ax 
shows  similarities  to  the  error  behavior  of  the  displacement  models  but  there  is  no  uniformity  in  the 
error  behavior.  For  example,  for  either  constant  Aror  Ax  the  error  values  reach  a  minimum  but  as 
time  progresses  these  minimum  error  values  do  not  follow  a  regular  pattern  as  for  the  displacement 
models.  Furthermore,  the  difference  in  error  between  the  CFET  and  the  CFEM  model  is  much  larger 
than  the  difference  between  the  DFET  and  DFEM  models. 

Another  observation  is  the  difference  in  error  between  the  displacement  and  conventional  models. 
In  all  six  figures  the  error  values  for  the  displacement  models  are  two  to  three  times  smaller  than  the 
corresponding  values  of  the  conventional  models.  For  example,  in  Figs.  2a  and  2b  at  small  values  of  Ar 
the  error  is  about  the  same,  but  as-Ar  increases  (Ar  >  0.01,  .V/,  >  .2)  the  error  for  the  conventional 
models.  Fig.  2b,  becomes  much  larger.  This  is  due  mainly  to  the  error  involved  in  applying  the  boun¬ 
dary  conditions  (x  =  L),  and  also  due  to  the  absence  of  the  boundary  forces.  In  the  case  of  the  con¬ 
ventional  models  the  boundary  forces  are  zero  for  prescribed  temperature  at  the  boundary.  Further¬ 
more,  as  was  mentioned  earlier  in  this  study,  the  presence  of  space  derivatives  in  the  complementary 
formulation  is  the  main  drawback  of  the  conventional  models  and  they  contribute  to  the  larger  error  of 
these  models.  In  the  following  discussion  certain  examples  will  be  used  to  demonstrate  this  defficiencv. 

The  previously  made  observations  on  the  patterns  of  the  error  behavior  suggest  that  optimum 
combinations  of  A  rand  Ax  can  be  defined  for  error  control  of  the  numerical  solution.  These  combina¬ 
tions  can  be  based  on  optimum  values  of  the  ratio  Xfx  =  Ar/Ax  or  the  ratio  .V/:  =  Ar/Ax2.  Typically 
for  a  second  order  equation,  such  as  the  transport  equation,  the  ratio  A/2  is  the  parameter  for  error  and 
stability  control.  The  value  of  this  parameter,  when  it  is  kept  within  certain  limits,  guarantees  accuracy 
and  stability  of  the  numerical  solution  [8], 
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A  careful  examination  of  the  results  for  the  error,  shown  in  Figs.  la-6a,  and  the  values  of  .V/; 
from  Tables  4  and  5  which  correspond  to  minimum  error  shows  that  there  is  an  inconsistency.  For 
example,  at  time  i  =*  0.5  Table  4  and  Fig.  2a  show  minimum  error  at  V/2  -  5.00,  (At  -  0.0125),  but 
from  Table  5  and  Fig.  5a  the  minimum  error  is  at  V/2  =  12.5  (,V£  =»  50).  What  appears  to  be  an 
inconsistency  in  the  behavior  of  the  displacement  models,  with  respect  to  the  parameter  A/2,  leads  to 
the  observation  that  Xl2  might  not  be  the  right  parameter  for  error  control  of  the  numerical  solution. 
In  contrast  an  examination  of  the  error  with  respect  to  A/,  shows  that  at  time  t  =  0.5  the  minimum 
values  of  the  error.  Figs.  2a  and  5a,  correspond  to  the  same  value  of  the  parameter  A/,.  In  fact  for 
each  value  of  the  parameter  \IX  the  corresponding  error  values  are  about  the  same  in  all  six  figures. 

In  order  to  study  the  dependence  of  the  numerical  solution  error  on  the  two  parameters  A/,  and 
.V/2  further  analysis  is  needed.  The  error  of  the  numerical  solution  now  is  investigated  for  some  con¬ 
stant  values  of  the  two  parameters.  Results  are  given  for  all  four  models  and  it  is  shown  that  uniform 
convergence  can  be  achieved  either  for  constant  M\  or  V/2.  Furthermore,  it  is  shown  that  uniform 
convergence  does  not  always  correspond  to  minimum  error. 

4.3  Convergence  of  the  Numerical  Solution 

Numerical  experiments  on  the  error  of  the  solutions  in  the  previous  section  show  that  nonuni¬ 
form  convergence  is  obtained  when  either  Ar  or  Ax  alone  is  the  dependent  variable.  It  was  also 
observed  that  minimum  error  could  be  achieved  for  some  optimum  combinations  of  the  two  variables 
A/ and  Ax.  It  is  a  well  known  fact  that  in  order  to  achieve  uniform  convergence  both  A  rand  Ax  should 
be  the  dependent  variables  with  the  condition  that  the  ratio  Ar/A.x  or  At/Ax2  is  kept  constant.  The 
results  for  the  average  error  £„  in  this  section  will  show  that  uniform  convergence  is  obtained  for  both 
the  displacement  and  conventional  models  when  M\  or  M2  is  kept  constant.  Some  representative 
values  of  Mx  and  A/2  are  given  in  Tables  8  and  9,  with  the  corresponding  values  of  Arand  Ax,  for  the 
four  models.  Numerical  results  are  presented  in  Figs.  7  through  15,  for  all  models  and  all  cases  as 
before,  and  for  two  values  of  each  of  the  parameters  M\  and  M2. 

Results  are  given  in  the  first  three  figures.  Figs.  7a-9a,  for  the  displacement  models  for  constant 

A/,  -  0.25  and  in  the  next  three.  Figs.  7b-9b,  results  are  given  for  M2  =  2.5.  The  uniform  conver- 
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Table  8 


V/, 

V/, 

A.v 

It 

2.50 

0.10000 

0.02500 

5.00 

0.0500 

0.01250 

0.25 

10.00 

0.0250 

0.00625 

12.50 

0.0200 

0.00500 

20.00 

0.0125 

0.00312 

25.00 

0.0100 

0.00250 

5.00 

0.1000 

0.05000 

0.50 

10.00 

0.0500 

0.02500 

20.00 

0.0250 

0.01250 

40.00 

0.0125 

0.00625 

50.00 

0.0100 

0.00500 

Table  9 


A/, 

V/2 

lx 

It 

0.250 

0.1000 

0.02500 

0.125 

0.0500 

0.00625 

0.063 

2.5 

0.0250 

0.00156 

0.050 

0.0200 

0.00100 

0.031 

0.0125 

0.00039 

0.025 

0.0100 

0.00025 

0.500 

0.1000 

0.05000 

0.250 

0.0500 

0.01250 

0.125 

5.0 

0.0250 

0.00312 

0.100 

0.0200 

0.00200 

0.063 

0.0125 

0.00078 

0.050 

0.0100 

0.00050 

gence  of  the  numerical  solution  is  obvious  for  either  of  the  parameters  M\  or  1 12.  Examining  each  pair 
of  figures  (Figs.  7a  and  b,  8a  and  b,  and  9a  and  b)  which  correspond  to  the  same  set  of  the  parameters 
V0  and  D0  one  can  observe  the  effect  of  /V/,  and  M2  in  the  error  of  the  solution  for  the  displacement 
models.  For  both  M\  and  M2  the  error  reaches  a  constant  value  but  these  constant  error  values  are 
smaller  for  the  case  of  constant  M\  than  for  the  case  of  constant  ,V/:.  This  becomes  more  evident  as 
the  value  of  D„  increases.  For  example,  in  Figs.  8a  and  8b,  for  each  value  of  lx  the  corresponding 
error  value  is  much  larger  for  constant  \f2  than  for  constant  A/,.  Furthermore,  the  corresponding 
value  of  It  for  each  lx,  Table  8,  is  less  for  the  case  of  M2  =  2.5  than  for  A/,  =  0.25.  This  behavior 
of  error  is  consistant  with  the  error  behavior  of  the  previous  section  where  it  was  shown  that  a  smaller 
value  of  It  for  a  particular  value  of  A*  does  not  always  mean  smaller  error. 
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Fig.  8  —  Average  error  for  the  displacement  models,  (  V0  -  1.0,  D0  -  1.0): 
(a)  A/,  -  0.25  and  (b)  Afj  -  2.5 
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Since  in  any  numerical  solution  we  are  concerned  not  only  with  the  convergence  but  also  with  the 
amount  of  error  involved  and  the  economy  with  which  the  solution  is  obtained,  it  is  obvious  that  in  the 


case  of  the  displacement  models  the  parameter  for  error  control  should  be  -V/,  and  not  A/:.  In  the  next 
six  figures.  Figs.  10-12  results  are  given  for  the  conventional  models  and  for  the  same  cases  as  above. 
Again  uniform  convergence  is  obtained  for  both  parameters  M\  and  Xf2  but  the  error  values  for  con¬ 
stant  A/,  are  now  larger  than  for  constant  .V/3.  This  is  the  opposite  from  what  was  the  case  for  the  dis¬ 
placement  models.  Furthermore,  the  difference  in  error  between  the  two  conventional  models,  CMET 
and  CFEM,  is  much  larger  for  A/,  =  0.25  than  for  A d2  =  2.5.  It  is  obvious  from  the  six  figures  that  for 
the  conventional  models  the  dominant  parameter  is  Xf2  and  not  A/,.  Also  for  the  larger  values  of  A t 
(.V/[  =  0.25)  the  numerical  solution  is  more  sensitive  to  the  type  of  time  integration  scheme  used. 

It  appears  that  for  the  displacement  models  the  correct  parameter  for  error  control  is  the  ratio  M\ 
and  for  the  conventional  models,  the  correct  parameter  is  the  ratio  \f2.  This  can  be  justified  first  by 
the  type  of  equation  involved  in  each  of  the  formulations  and  second  by  the  results  obtained  for  the 
average  error  of  the  numerical  solutions.  The  equation  for  which  numerical  solutions  are  sought 
through  the  displacement  models  Eq.  (8),  does  not  contain  second  order  space  derivatives  of  the  tran¬ 
sport  variable  and  only  the  ratio  Ar/A.v  appears  in  the  approximations  of  this  equation.  In  contrast  the 
equation  for  which  numerical  solutions  are  sought  through  the  conventional  models,  Eq.  (3),  does  con¬ 
tain  second  order  derivatives,  hence  the  ratio  A&'A.v2  appears  in  the  approximation  of  this  equation. 
The  results  presented  in  this  section  are  in  full  agreement  with  this  reasoning. 

In  the  last  part  of  this  section  two  more  sets  of  figures  are  given  for  different  values  of  the  param¬ 
eters  M ]  and  A/,.  These  sets  of  figures  together  with  Figs.  7a-9a  and  Figs.  10b- 12b  can  be  used  to 
evaluate  and  compare  the  displacement  and  conventional  models.  For  the  displacement  models.  Figs. 
7a-9a  and  Figs.  1 3a- 1 5a  the  error  increases  as  the  value  of  A /t  increases  from  0.25  to  0.5.  A  similar 
behavior  is  observed  for  the  conventional  models  as  M2  increases  from  2.5  to  5.0.  Comparing  now 
Figs.  7a  through  9a  to  Figs.  10b  through  12b  and  Figs.  13a  through  15a  to  Figs.  13b  through  15b  one 
can  observe  that  the  displacement  models  consistently  result  in  smaller  error  values  than  the  conven- 


tionai  ones,  especially  for  large  values  of  Ax  (small  NE).  Furthermore,  the  smaller  error  values  are 
obtained  by  the  displacement  models  with  a  larger  time  step.  For  a  given  value  of  Ax,  say  Ax  =  0.025, 
the  values  of  Ar  which  correspond  to  .V/,  -  0.25  (DFET  and  DFEM)  and  M2  —  2.5  (CFET  and  CFEM) 
are  Ar*  =  0.00625  and  A t2  =  0.00156  respectively.  This  can  be  easily  observed  from  the  two  Tables  8 
and  9  and  from  Figs.  7  through  15. 

The  following  general  observations  can  be  made  regarding  the  error  behavior  of  the  numerical 
models. 

•  Uniform  convergence  is  obtained  for  all  models  only  when  the  ratio  M\  or  the  ratio  A/2  is  kept  con¬ 
stant. 

•  Uniform  convergence  does  not  automatically  guarantee  minimum  error. 

•  For  the  DFET  and  DFEM  models  the  minimum  error  is  obtained  when  the  control  parameter  is 

Mx. 

•  For  the  CFET  and  CFEM  models  the  minimum  error  is  obtained  when  the  control  parameter  is 

.Vf2. 

•  Of  the  two  time  integration  schemes  used,  the  finite  element  in  time  is  consistantly  more  accurate 
in  all  cases. 

•  The  displacement  models  are  computationally  more  efficient  overall  since  they  can  achieve  smaller 
error  than  the  conventional  ones  by  using  larger  time  steps. 

•  For  all  of  the  models  the  error  decreases  significantly  up  to  A '£  =  40  and  an  increase  of  NE  beyond 
this  value  has  only  small  effect  on  the  numerical  solution. 
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Fig.  10  —  Average  error  for  the  conventional  models,  <  F0  -  1.0.  D0  -  0.001): 

(a)  -  0.25  and  (b)  M7  -  2.5 
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Fig.  12  —  Average  error  for  the  conventional  models.  ( F0  -  0.0.  D0  -  1.0): 
(a)  ,W,  -  0.25  and  (b)  ,Vf2  -  2.5 
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Fig.  14  —  Average  error  (  V0  -  1.0.  D0  -  1.0):  (a)  Displacement  models 
A/,  —  0.5  and  (b)  Conventional  models  M2  -  5.0 
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Numerical  experimentations  based  on  the  average  error  result  not  only  in  a  good  understanding  of 
the  overall  behavior  of  the  numerical  models  but  also  how  the  numerical  solutions  obtained  by  these 
models  are  affected  by  the  parameters  .V/,  or  Xl2.  From  the  results  presented  in  the  last  two  sections 
and  from  previous  studies  some  optimum  values  for  A/,  and  M2  can  be  determined.  These  optimum 
values  should  correspond  to  small  error  but  at  the  same  time  should  hold  the  efficiency  of  the  computa¬ 
tional  mode!  at  a  high  level.  The  range  of  values  for  .V/,  is  between  0.1  and  0.3  and  for  >/2  is  between 
2.0  and  4.00.  The  exact  choice  of  value  sometimes  depends  on  the  conditions  of  the  specific  problems 
and  one  might  have  to  choose  a  smaller  value  for  M\  or  \f2. 

The  experience  gained  from  studying  the  error  behavior  can  be  used  for  obtaining  solutions  to 
specific  problems.  From  such  solutions  one  can  evaluate  the  efficiency  of  a  computational  model  in  a 
more  detailed  fashion,  and  with  respect  to  other  parameters  such  as  boundary  conditions.  In  the  next 
section  several  problems  will  be  presented  for  that  purpose. 

4.4  Numerical  Examples 

The  computational  models  developed  in  this  study  will  be  applied  to  solve  the  general  transport 
equation  for  a  number  of  physical  problems.  Numerical  solutions  of  the  transport  equation,  with  the 
associated  boundary  conditions  are  discussed  in  this  section  and  results  obtained  from  the  displacement 
(DFET)  and  conventional  (CFET)  models  are  compared  to  available  analytical  solutions.  Results  for 
all  case  studies  were  obtained  by  using  XE  =  50  (A.v  »  0.02)  and  .V/,  =  0.25  for  the  DFET  model  and 
Xf2  =*  2.5  for  the  CFET  one.  The  corresponding  time  steps  are  indicated  on  the  figures  for  each  case 
study. 

4.4.1  Advection  Equation 

This  is  the  simplest  form  of  the  general  transport  equation  representing  a  number  of  physical 
processes  such  as  mass  conservation,  convection  of  heat  and  free  surface  elevation  in  flow  problems. 


50 


f>(x,o )  =  0.0  0.0  <  x  ^  1.0 

p(o.t)  =1.0  t  >  0.0. 


(95) 


Numerical  solutions  are  obtained  for  Eq.  (93),  displacement  models,  and  for  Eq.  (92).  conventional 
models.  For  the  displacement  model  the  density  p  is  evaluated  through  the  constraint  given  by  Eq. 
(94)  and  the  same  equation  is  used-  to  express  the  initial  and  boundary  condition  in  terms  of  the  dis¬ 
placement  field. 


The  first  two  figures  (16a  and  b)  show  results  obtained  by  the  two  models  for  the  classical  prob¬ 
lems  of  a  discontinuity  advected  through  space  at  four  different  time  levels.  Both  models  show  the 
same  typical  overshoots  and  the  overall  error  is  about  the  same  for  both  models.  No  attempt  was  made 
to  suppress  these  oscillations  since  a  realistic  evaluation  of  the  models  is  sought.  In  the  second  two 
figures  (17a  and  b)  results  are  given  for  a  different  case  of  the  advection  equation.  A  point  source  term 
is  now  included  in  the  equation,  with  constant  rate  of  production  equal  to  one.  The  source  is  placed  at 
x  —  0.5  and  the  mass  produced  is  advected  to  the  right  with  velocity  V0  -  1.0.  For  the  CFET  model 
some  oscillations  exist  in  the  region  (x  <  0.5)  before  the  source.  There  oscillations  are  much  smaller 
for  the  DFET  model.  The  oscillatory  behavior  of  the  CFET  model  is  due  to  the  discontinuity  in  the 
density  field  imposed  by  the  presence  of  the  source  term  at  x  =  0.5.  At  the  same  point,  the  displace¬ 
ment  field  is  continuous  and  has  a  slope  of  one  according  to  Eq.  (9).  The  absence  of  discontinuities  in 
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the  displacement  field  not  only  reduces  the  oscillations  about  x  =  0.5  but  also  results  in  a  more  stable 
solution,  as  one  can  observe  from  the  values  of  the  density,  close  to  the  left  hand  boundary  (x  =  0). 
Some  of  the  solutions  in  the  next  examples  show  similar  behavior  and  by  a  closer  look  at  the  displace¬ 
ment  model  one  can  observe  that  boundary  conditions  are  better  handled  when  they  are  imposed 
through  the  displacement  formulation. 

Although  the  advection  equation  is  the  simplest  form  of  partial  differential  equation,  its  numerical 
solution  is  challenging  and  a  large  number  of  publications  exist  on  other  types  of  initial  or  boundary 
conditions  as  ’-'ell  as  different  approaches.  In  the  present  cases  the  finite  element  method,  in  its  sim¬ 
plest  form,  ,*as  applied  for  the  two  derived  formulations.  Since  the  main  objective  here  is  to  evaluate 
the  conventional  and  displacement  formulations,  we  will  continue  to  use  linear  finite  element  approxi¬ 
mations  for  all  cases  of  the  transport  equation. 


4.4.2  Convection  —  Diffusion 


The  inclusion  of  a  second  order  space  derivative  in  the  advection  equation  results  in  the 
convection-diffusion  equation.  This’  is  the  governing  equation  for  a  number  of  physical  processes  such 
as  convective  heat  transfer  and  energy  balance,  convective  transport  of  species,  pollutants,  and  momen¬ 
tum  transfer.  The  basic  convection-diffusion  equation  is  given  by 


for  a  constant  velocity  field  and  diffusivity. 


Here  the  transport  variable  is  represented  by  the  dimensionless  temperature  Q(x,r )  which  relates 
to  the  transport  (heat)  displacement  in  the  same  fashion  as  for  the  density 
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dH(x.r) 

dx 


(98) 


If  a  source  term  and  a  heat  transfer  term  are  included  in  the  convection-diffusion  equation,  then  Eq. 


(97)  is  written  as 


i®.  +  v  -z> 

d/  0  0*  0 


d20 

d-v2 


S0S  -  Ka& 


where  S0  is  the  coefficient  associated  with  the  heat  source  term  S  ( S0 


(99) 

0  or  1)  and  Kn  is  the  heat 


transfer  coefficient. 


In  terms  of  the  heat  displacement  field  H(x,t )  the  convection-diffusion  equation  is  written  as 

S„  h  -  KnH  (100) 


dH  ,  ML  -  n  OIK 
dr  0  dx  0  dx2 


Initial  and  boundary  conditions  associated  with  this  equation  are  numerous.  Here  some  typical  boun¬ 
dary  conditions  are  considered  and  numerical  solutions  will  be  obtained  for  both  computational  models. 
In  Table  10  the  initial  and  boundary  conditions  are  given  for  each  of  the  problems  considered. 


Table  10 


•So 

K0 

Initial  Conditions 

Boundary  Conditions  L  —  <* 

Figures 

1 

0.0 

0(1,0)  =  l.O.0(£.O)  -  0.0 

18 

2 

0.0 

0UO)  -  0.0 

0(1,0)  =  sin  a»r,0(Z.,O)  =  0.0 

19  ! 

3 

1.0 

1.0 

1.0 

0.0 

0(.v,O)  =  0.0 

0(  O.r)  =00,0(1,0)  -  0.0 

a 

4 

1.0 

1.0 

^ -  (0,/)  =O.O.0(1,O)  =  0.0 
dx 

5 

1.0 

1.0 

Q(-L.t)  =  0.0, 0(1,/)  =  0.0 

22 

6 

1.0 

1.0 

1.0 

1.0 

0(.v,O)  =  0.0 

0(0./)  =  0.0. 0(1,0)  =  0.0 

LJ1 _ 1 

In  the  first  problem  the  right  hand  side  boundary  condition  .(x  =  1)  is  assumed  to  be  applied  at  infinity. 
This  condition  will  be  assumed  in  all  problems  considered  in  the  following  examples.  For  the  heat 
source  term  the  value  of  S„  is  1.0  at  x  =  0.5  and  0.0  elsewhere  with  li(O.S.r)  =  1.0  This  represents  a 
point  heat  source  at  .v  =  0.5  with  unit  heat  generation  per  unit  time.  Numerical  results  for  the  prob¬ 
lems  described  in  Table  10  are  given  in  Figs.  18  through  23. 

The  first  problem  is  the  typical  case  for  the  convection-diffusion  equation  for  a  constant  tempera¬ 
ture  applied  on  the  left  hand  side  boundary.  The  physical  boundary  to  the  right  is  at  infinity  and  the 
computational  boundary  is  at  x  =  1.0.  Since  the  boundary  condition  0(1./)  *  0.0,  is  at  infinity,  then  it 
is  physically  incorrect  to  apply  this  boundary  condition  at  .v  =  1.0.  The  temperature  at  x  -  1.0  is  zero 
for  small  times  (t  <  1.0)  but  it  is  not  zero  for  larger  times.  Therefore  no  boundary  condition  is 
imposed  at  x  =  1.0  and  the  computational  model  should  be  able  to  simulate  an  outflow  boundary  at 


that  point. 
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In  Fig.  18a  results  for  the  DFET  model  are  given  for  four  time  intervals  and  the  numerical  solu¬ 
tion  is  compared  to  the  analytical  one  (solid  lines).  The  agreement  between  the  two  is  good  and  the 
largest  error  occurs  close  to  the  right  hand  side  boundary.  Similar  results  are  given  in  Fig.  18b  for  the 
CFET  model  for  the  temperature  distribution.  The  error  of  the  CFET  model  is  larger  than  the  DFET 
one  throughout  the  computational  domain.  The  larger  error  occurs  mainly  because  of  the  error  intro¬ 
duced  by  the  boundary  conditions.  The  boundary  forces  resulting  from  the  variational  formulation  are 
given  by 


for  the  CFET  model,  and 


a©  a©  , 
s7  si* 


(101) 


(?,=*/  n©  * 
Js  a?, 

for  the  DFET  model. 


(102) 


If  a  boundary  condition  is  described  for  the  temperature  at  x  =*  0.0  (©(0,r)  —  1.0)  then 

se  -  0  —  |^-  -  0 
3  <h 

which  means  that  the  boundary  force  at  x  =  0.0  is  zero  for  the  CFET  model.  This  is  not  the  case  for 
the  DFET  model  since  8H,  is  not  zero.  At  the  other  boundary  <x  —  1.0)  the  same  is  true  for  both 
models  if  the  temperature  is  described  at  that  boundary.  As  was  mentioned  before  no  boundary  condi¬ 
tion  is  imposed  at  x  =  1.0,  which  means  that  the  boundary  forces  are  not  zero  and  should  be  included 
in  the  calculation  for  both  models.  By  retaining  the  boundary  forces  the  results  for  the  DFET  model. 
Fig.  18a,  show  a  good  agreement  with  the  analytical  solution  and  the  boundary  force  accounts  for  the 
neglected  portion  of  the  physical  domain  (x  >  1).  The  inclusion  of  the  boundary  force  in  the  compu¬ 
tations  for  the  DFET  model  is  also  in  agreement  with  the  theory  since  neither  ©d.r)  nor  5 H,  is  zero. 

Although  theory  and  results  are  in  agreement  for  the  displacement  formulation  the  same  is  not 
true  for  the  conventional  one.  If  one  retains  the  boundary  force  in  the  compulations,  the  results,  given 
in  Fig.  18c,  show  very  large  error  close  to  the  boundary  (x  -  1.0)  but  if  the  boundary  force  is  set  equal 
to  zero,  even  though  0  or  8©  are  not  zero,  then  the  results,  given  in  Fig.  18b,  show  much  less  error 
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Fig.  19  —  Temperature  distribution  (  F0  -  1.0,  D0  -  1.0.  Problem  2): 

(a)  Displacement  model  and  (b)  Conventional  model 
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and  the  temperature  distribution  is  similar  to  the  DFET  model.  This  discrepancy  in  the  results 
obtained  by  the  CFET  model  exists  in  all  cases  of  solving  problems  in  the  semi-infinite  space.  There  is 
a  deficiency  of  the  conventional  formulation  in  terms  of  properly  simulating  outflow  conditions.  For  all 
these  types  of  problems  even  if  the  boundary  conditions  are  not  explicitly  imposed,  the  boundary  forces 
should  be  zero  for  the  CFET  model.  This  is  in  agreement  with  the  Galerkin  approach  and  the  term 
"consistent  model"  will  be  used  for  the  CFET  model  with  no  boundary  forces.  For  the  second  problem, 
where  a  sinusoidal  wave  (to  =  tt/,2  )  is  applied  at  the  boundary  (x  =  0),  results  are  presented  in  Figs. 
19a  and  19b  for  each  of  the  two  models.  A  comparison  of  these  two  figures  shows  a  good  arguement 
between  the  two  models  but  it  should  be  noted  that  the  results  for  the  CFET  model  were  obtained  by 
using  a  time  step  five  times  smaller  than  the  one  used  for  the  DFET  model.  This  is  the  case  for  all  the 
results  obtained  in  this  study,  where  in  order  to  achieve  results  of  about  the  same  accuracy  the  CFET 
model  requires  a  time  step  five  times  smaller. 

The  next  example,  problem  3,  is  the  case  when  a  point  heat  source  is  present  at  x  =  0.5.  Results 
for  the  DFET  model  are  given  in  Fig.  20a  and  for  the  CFET  model  in  Figs.  20b  and  20c.  In  the  Jasi 
two  figures  results  are  given  for  cases  of  omission  and  inclusion  of  the  boundary  force  (x  =—  1.0) 
respectively  as  previously  for  problem  1.  As  one  can  observe  from  Fig.  20b  (no  boundary  force)  the 
results  for  the  consistent  model  are  in  good  agreement  with  the  ones  of  the  DFET  model.  Fig.  20a. 
But  when  the  boundary  force  is  included.  Fig.  20c,  the  temperature  distribution  results  in  much  larger 
values  close  to  the  boundary  (x  =*  1.0).  This  behavior  of  the  CFET  model  is  the  same  as  previously  in 
problem  1.  For  the  rest  of  the  problems  in  this  section  the  results  for  the  CFET  model  are  with  no 
boundary  force.  In  problem  4  a  point  heat  source  is  present,  as  in  the  previous  problem,  and  the  boun¬ 
dary  conditions  are  of  no  flow  of  heat  at  the  left  boundary  and  zero  temperature  at  infinity.  The  results 
for  the  two  models  are  presented  in  Fig.  21a  and  21b.  The  temperature  distribution  shows  good  agree¬ 
ment  at  the  left  boundary  but  again  at  the  outflow  boundary  the  temperature  is  lower  for  the  CFET 
model.  This  has  been  the  case  for  all  three  previous  examples.  If  one  now  considers  a  infinite  domain 
and  that  the  temperature  is  zero  at  -L  and  +  L  (L  —  <»),  the  heat  generated  by  a  point  heat  source 
will  be  convected  to  the  right  and  diffused  in  both  directions  according  to  the  governing  equation. 
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Fig.  20  —  Temperature  distribution  ( yg  -  1.0,  D0  -  1.0,  Problem  3): 
(a)  Displacement  model,  (b)  Conventional  model  and  (c)  Conventional 
model  with  boundary  forces 
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Fig.  21  —  Temperature  distribution  ( F0  —  1.0.  Do”  ' -0.  Problem  4): 
(a)  Displacement  model  and  lb)  Conventional  model 


Numerical  results  for  the  solution  of  this  problem  are  given  in  Fig.  22a  and  22b  for  the  two  models. 
The  temperature  distribution  obtained  from  the  CFET  model  is  almost  twice  as  high  than  the  one 
obtained  in  the  DFET  model  at  the  left  side  of  the  computational  boundary  and  it  is  lower  to  the  right 
of  the  heat  source.  Again  this  is  due  to  the  way  that  outflow  boundary  conditions  are  simulated  by  the 
CFET  model.  As  a  last  example  for  the  convection  and  diffusion  of  heat  a  case  is  presented,  problem 
5,  where  heat  is  lost  throughout  the  space  ((0  <  x  <  °°))  due  to  the  heat  transfer  term  in  the  govern¬ 
ing  equation.  Results  for  this  case  are  given  in  Fig.  23a  and  23b  for  the  two  models.  These  two  figures 
should  be  compared  with  Figs.  20a  and  20b  where  results  are  given  for  the  same  problem  but  with  no 
heat  loss.  As  it  is  expected  the  temperature  distribution  for  problem  5  is  lower  than  the  one  for  prob¬ 
lem  3  throughout  the  computational  domain. 

For  all  case  studies  considered  here  both  computational  models  give  good  results  for  problems 
where  the  boundary  conditions  are  well  defined.  On  the  other  hand  when  the  outflow  conditions  exist, 
the  CFET  model  yields  results  at  that  boundary  which  are  always  lower  to  the  corresponding  results  of 
the  DFET  model.  Furthermore,  for  such  cases  where  analytical  solutions  are  available,  a  comparison 
with  the  numerical  solutions  shows-  that  the  CFET  model  yields  always  a  lower  distribution  than  the 
DFET  model. 

The  failure  of  the  CFET  model  to  simulate  correctly  outflow  conditions  is  related  mainly  to  the 
fact  that  the  conventional  model  cannot  account  for  the  presence  of  the  physical  space  beyond  the  com¬ 
putational  one.  This  deficiency  results  in  more  heat  flowing  out  at  the  end  of  the  computational  boun¬ 
dary,  thus  lower  temperature  distributions.  In  the  next  section  similar  behavior  is  observed  for  the 
diffusion  of  heat  in  the  half  space. 

4.4.3  Diffusion  Equation 


In  the  case  of  a  solid  medium  or  a  fluid  not  in  motion  the  convective  term  of  the  transport  equa¬ 
tion  is  zero  and  the  governing  equation  is  given  by 

a2© 


— —  D 
dt  "  Bx2 
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(103) 
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Fig.  23  —  Temperature  distribution  (  P0  -  1.0,  £>0  -  1.0,  Problem  6): 
(a)  Displacement  model  and  (b)  Conventional  model 
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with  0(.v.;)  the  dimensionless  temperature  which  is  related  to  the  heat  displacement  Hix.t)  by  Eq. 


Problem 

A, 

A 

A'„ 

Initial 

Conditions 

Boundary 

Conditions 

Figure 

1 

1.0 

0.0 

0.0 

0<x,O)  =  0.0 

0(0,;)  =  1.0 

0(1.;)  -  0.0 

24 

2 

1.0 

0.0 

0.0 

0(x,O)  =  0.0 

0(0,;)  —  sin(tur) 
©(£,;)  =  0.0 

25 

3 

1.0 

1.0 

0.0 

0(.v.  0)  =  0.0 

0(0.;)  =  0.0 

0(1.;)  =  0.0 

26  1 

4 

1.0 

1.0 

0.0 

0(.v.  0)  =  0.0 

0(— I.;)  =  0.0 
©(A;)  =  0.0 

27 

5 

1.0 

0.0 

0.0 

0(.v,  0)  =  0.0 

0(0.;)  =  1.0 

0(1.;)  -  0.0 

28 

6 

1.0 

0.0 

0.0 

0(.v.  0)  =  0.0 

0(0,;)  -  1.0 

^ -  (1.;)  =  0.0 

8.v 

29 

7 

1.0 

0.0 

0.0 

0(x,  0)  =  0.0 

0(0,;)  =  1.0 

d.t)  =  0(U) 

ox 

30 

8 

1.0 

1.0 

0.0 

0(xO)  =  0.0 

0(0,;)  =  0.0 

0(1.;)  =  0.0 

31 

9 

1.0 

1.0 

0.0 

0(.v,O)  =  0.0 

-00  1 
self.;)  „00 

dx 

32 

10 

1.0 

1.0 

1.0 

0(.v.O)  =  0.0 

-mrnirr^r 

»<*'(•'»  -00 
dx 

33 

The  initial  and  boundary  conditions  associated  with  Eq.  (103),  for  the  case  studies  considered 
here,  are  given  in  Table  11.  The  first  four  problems  are  for  the  diffusion  of  heat  in  the  semi-infinite 
space  (£—«>)  and  the  last  six  for  problems  of  finite  space  (0.0  <  x  <  1.0)  .  The  first  problem  is  the 
classical  case  of  constant  temperature  at  .v  =  0.0  and  zero  temperature  at  infinity.  Results  are  given  in 
Figs.  24a  and  24b,  for  the  two  models,  and  they  are  compared  to  the  analytical  solution  (solid  line).  As 
was  previously  observed,  for  the  case  of  convection-diffusion,  the  temperature  distribution  for  the 
CFET  model  is  lower  than  the  DFET  model.  The  error  in  the  results  of  the  CFET  model  is  larger  for 
the  diffusion  case  than  for  the  convection  diffusion  one.  This  should  be  expected  since  the  diffusion 
term  is  the  dominent  term  in  the  present  case,  thus  the  boundary  force  plays  a  more  important  role  in 
simulating  outflow  conditions.  Similar  observations  can  be  made  for  the  third  and  fourth  problems. 
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Fig.  24  —  Temperature  distribution  f  V0  —  0.0,  Dv  —  1.0,  Problem  1): 
(a)  Displacement  model  and  <b)  Conventional  model 
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Fig.  25  —  Temperature  distribution  (  V0  -  0.0,  D0  -  1.0.  Problem  2): 
(a)  Displacement  model  and  (b)  Conventional  model 
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Figs.  26  and  Figs.  27.  For  both  problems  a  point  heat  source  is  located  at  x  —  0.5  and  the  boundary 
conditions  are  given  in  Table  11.  For  problem  3  a  comparison  of  results  between  the  DFET  model. 
Fig.  26a  and  the  consistent  CFET  model  (no  boundary  force  at  x  =  1.0),  Fig.  26b,  shows  a  much  lower 
temperature  distribution  to  the  right  of  the  heat  source.  The  opposite  is  true  for  the  results  of  the 
CFET  model  of  Fig.  26c. 

In  the  next  problem  where  the  computational  domain  represents  an  infinite  physical  domain  the 
temperature  distribution  should  be  symmetric  with  respect  to  the  point  heat  source.  This  is  true  for  the 
results  of  the  DFET  model.  Fig.  27a,  and  the  results  of  the  CFET  model  when  the  boundary  forces  are 
included.  Fig.  27c  but  it  is  not  true  for  the  results  obtained  by  the  consistant  CFET  model.  Throughout 
the  results  presented  in  the  previous  and  present  sections  an  important  difference  exists  between  the 
two  models.  The  implementation  of  the  DFET  model  to  solve  problems  with  unbounded  domains 
results  in  computational  models  where  no  inconsistancies  exist  between  the  physical  and  mathematical 
aspects  of  the  problem.  Although  this  is  the  case  for  the  DFET  model,  it  is  not  true  for  the  CFET 
model  where  the  absence  of  the  boundary  forces  result  in  a  better  numerical  solution  but  the  physical 
and  mathematical  aspects  of  the  problem  contradict  each  other. 

In  the  last  part  of  this  section  the  diffusion  of  heat  is  considered  for  a  finite  domain,  a  slab  of 
thickness  one,  and  the  boundary  condition  for  the  different  cases  are  given  in  Table  11.  In  all  exam¬ 
ples,  problems  5-10,  the  numerical  results  of  the  two  models  are  of  about  the  same  accuracy.  For  the 
finite  domain  problems,  the  boundary  conditions  are  applied  on  both  boundaries.  The  boundary  forces 
are  zero  for  the  CFFT  model  (all  examples)  but  for  the  DFET  model  the  boundary  forces  are  zero  only 
if  the  specified  boundary  temperature  is  zero.  Otherwise  the  values  of  the  boundary  forces  depend  on 
the  values  of  the  temperature  at  each  boundary. 

The  numerical  results  of  problems  5  through  10  are  presented  in  Figs.  28  through  33  for  the  two 
models.  Comparing  the  results  obtained  by  the  two  models  one  can  easily  observe  that  the  numerical 
solutions  are  almost  identical.  The  only  difference  is  that  the  results  of  the  CFET  model  were  obtained 
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Fig.  26  —  Temperature  distribution  (V0  —  0.0.  O0  —  1.0,  Problem  3): 
(a)  Displacement  model,  (b)  Conventional  model  and  (c)  Conventional 
model  with  boundary  forces 
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Fig.  27  —  Temperature  distribution  ( yg  -  0.0.  Do«1.0,  Problem  4): 
(a)  Displacement  model,  (b)  Conventional  model  and  (c)  Conventional 
model  with  boundary  forces 
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Fig.  31  —  Temperature  distribution  (K#—  0.0,  D0  —  1.0,  Problem  8): 

(a)  Displacement  model  and  (b)  Conventional  model 
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Fig.  32  —  Temperature  distribution  ( -  0.0,  Dg  -  1.0,  Problem  9): 

(a)  Displacement  model  and  (b)  Conventional  model 
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Fig.  33  —  Temperature  distribution  (  F0  —  0.0,  D0  —  1.0,  Problem  10): 
(a)  Displacement  model  and  (b)  Conventional  model 


by  using  a  time  step  size  five  times  smaller  than  the  one  used  for  the  DFET  model.  Therefore  from 
the  computational  point  of  view  the  displacement  formulation  is  more  efficient  since  it  requires  less 
time  to  solve  a  problem.  If  the  problem  to  be  solved  does  not  have  outflow  boundary  conditions,  both 
models  produce  numerical  solutions  of  the  same  accuracy. 

Throughout  the  examples  presented  here  one  can  observe  the  displacement  formulation,  the 
DFET  model,  is  more  consistent  than  the  conventional  one  with  respect  to  handling  different  physical 
problems.  This  is  justified  not  only  from  the  numerical  solutions  obtained  but  also  from  the  formula¬ 
tion  on  which  the  displacement  model  is  based. 

SUMMARY  AND  CONCLUSIONS 

A  variational  formulation  for  the  transport  equation  has  been  presented  in  this  study,  and,  based 
on  this  formulation,  a  generalized  approach  to  approximate  solutions  has  been  developed.  As  a  special 
application  to  approximate  solutions  four  computational  models  were  derived  for  obtaining  numerical 
solutions  to  different  types  of  transport  problems.  The  models  are  based  on  first-order  finite  element 
approximations  for  both  the  spatial  arid  time  coordinates.  Two  of  the  models  are  derived  from  the  fun¬ 
damental  variational  analysis  in  terms  of  the  transport  displacement  and  the  other  two  are  derived  from 
the  complementary  variational  analysis  in  terms  of  primitive  transport  variable.  Although  only  linear 
and  one-dimensional  models  were  considered  here,  higher  order  and  multi-dimensional  models  can  be 
easily  obtained  since  the  derived  formulation  is  not  restricted  to  one  type  of  approximation  or  one 
dimension. 

An  evaluation  of  the  four  computational  models  with  respect  to  accuracy  and  efficiency  was  car¬ 
ried  out,  first,  by  numerical  experimentation  on  the  error  behavior  and,  second,  by  obtaining  numerical 
solutions  to  specific  boundary  value  problems. 

The  investigation  of  the  error  of  the  numerical  solution  shows  that  the  displacement  models  not 
only  are  more  consistent  in  the  error  behavior  than  the  conventional  ones,  but  also  that  the  overall 
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error  of  the  numerical  solutions,  obtained  by  the  displacement  models,  is  always  smaller  than  the  error 
of  the  solutions  obtained  by  the  conventional  ones.  The  larger  error  of  the  conventional  models  is  par¬ 
tially  because  they  are  derived  from  the  complementary  variational  principle,  which  according  to  the 
theory  is  less  accurate  than  the  corresponding  fundamental  variational  principle.  Numerical  solutions  to 
heat  transport  problems  obtained  here  also  justify  this  argument.  Furthermore,  these  solutions  showed 
the  inability  of  the  conventional  models  to  correctly  simulate  outflow  boundary  conditions. 

Although  the  two  types  of  computational  models  have  some  common  characteristics  in  behavior, 
the  conventional  ones  are  less  efficient  from  both  accuracy  and  computational  effort  points  of  view. 
Displacement  models  have  been  used  widely  in  structures,  but  as  it  has  been  shown  these  models  can 
be  extended  to  other  types  of  problems  as  well.  Furthermore,  the  fundamental  form  of  the  variational 
analysis  can  be  derived,  based  on  concepts  of  classical  mechanics  presented  here,  for  many  governing 
equations. 

In  conclusion,  the  generalized  variational  analysis  and  the  unified  approach  in  obtaining  approxi¬ 
mate  solutions  should  be  emphasized  as  well  as  the  applicability  of  the  analysis  to  other  physical  prob¬ 
lems.  It  should  be  also  noted  that  the  finite  element  method,  which  was  presented  in  this  study,  is  only 
one  restricted  application  of  the  generalized  formulation. 
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